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ABSTRACT 


One dimensional flow between two fixed parallel walls composed 
of the same substance but at different temperatures and spaced a 
distance l apart is considered. The hot plate is the evaporating 
surface (source) and the cold plate is the condensing surface (sink). 
The vapor between the two plates is assumed to be a monatomic gas 
consisting of Maxwell molecules. Lees' moment method is used to 
obtain a set of six non-linear equations whose solution, subject to the 
boundary conditions of this problem, is possibly valid from free 
molecular to continuum conditions. 

Both the non-linear equations and a linearized approximation to 
them are solved. 

The non-linear problem required the solution of six simultan- 
eous and ordinary non-linear differential equations with three bound- 
ary conditions given at each wall. An iterative numerical procedure 
was used to match these boundary conditions. For the continuum 
limit (Reynolds number large), the vapor leaving the hot plate was 
found to accelerate rapidly to an equilibrium velocity. In the vicinity 
of the cold wall, the vapor first decelerated, then experienced a 
slight terminal acceleration. In the rarefied limit (Reynolds number 
very small), the vapor velocity was found to be essentially constant 

iii 


across the flow field. 



The linearized problem in closed form under the assumption 
of a small mean velocity is solved. Large and small Knudsen 
numbers are examined. In both the rarefied and continuum limits, 
the mean velocity was found to be constant across the flow field. 

For given emission temperatures and density ratios at the two 
surfaces, the mean speed between the plates varied with Reynolds 
number because of the effects of molecular collisions. 

The evaporation coefficient is defined here as the ratio of the 
actual mass flux to the difference between the Knudsen effluxes from 
the two surfaces. Its value is nearly one for the range of Knudsen 
numbers considered in both the linear and non-linear problems. 
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NOMENCLATURE 


B = mass constant 

= momentum constant 
B^ = energy constant 
c. = vector particle velocity 
C. = relative particle velocity 
C = mean molecular speed 
d = distance between two plates 

E X ,E 2 = 1 + erf j^u 1 /v / 2RT 1 ] . 1 + erf [-u^/TrT^ ] 

f,f^ = velocity distribution functions for "probe" and colliding 
particles, respectively 

f ,f •= components of two stream Maxwellian 

X £ 

F. = external force acting on a single particle 

h = specific enthalpy 

H = total enthalpy 

k = Boltzmann constant 

Kn = Knuds en number (X/d) 

m = molecular mass 

M - Mach number (u/a) 

n -- number of molecules per unit volume 

n^.n^ = number density functions in two stream Maxwellian 
N, = small perturbation to n 

1 j E i j u 

viii 




N_ = N r N 2 


N+ *VN 2 


p = pressure (p = pRT) 
q = heat flux vector 

Q = arbitrary function of particle velocity 
R = gas constant 
Re - Reynolds number 
S = speed ratio (u/ / 2RT ) 
t = time 

t = small perturbations to T 
^ ^ 1 1 2 

fc - = V fc 2 
V 4^2 

T = temperature 

T 1’ T 2 = temperature functions in two stream Maxwellian 
u = mean velocity in x direction 

U 1’ U 2 = velocit y functions in two stream Maxwellian 


U = /"RTj. or /RT 

U ,U = small perturbations to u ,u 
i ^ 12 

x = coordinate in x direction 
X,,X 2 = exp[- u ^ 2 /2RT 12 ] 

O' = the evaporation coefficient 


B l,2 = 2RT ll2 
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e = 

X = 

H = 
v = 

V l,2 
P = 

a. . = 
ij 


small parameter 
mean free path 
viscosity 

kinematic viscosity 

= u 12 //2F 

gas density 
stress tensor 


T. . = o . .+ p 6 . . 

1J 1J 1J 

Subscripts 

I indicates hot wall 


II indicates cold wall 



Chapter 1 


INTRODUCTION 

When the state of the vapor at a gas -liquid interface is at the 
saturated temperature and pressure corresponding to surface 
temperature, an equilibrium situation exists. There is an exact 
balance between the two molecular processes of evaporation and 
condensation. If the vapor arid liquid surfaces are not in equilibrium, 
then, a net flux of molecules either condenses on or evaporates from 
the surface. 

In recent years many publications have appeared dealing with 
evaporation from and condensation on a solid or liquid surface. 
Evaporation either into a near vacuum or into a gas where there are 
small deviations from equilibrium at the liquid-vapor interface have 
been the two areas most thoroughly studied. Little effort has been 
directed towards investigating the evaporation and condensation 
process over a wide pressure ratio range. 

Prior to the space program, interest in low pressure evapora- 
tion and condensation grew from the study of thin film deposition, 
molecular distillation, and other aspects of vacuum technology. In 
the last dec.aae the field has expanded to include the study of numer- 
ous materials which may evaporate or sublime in a space environ- 
ment. Such processes produced both forces and heat fluxes at the 

1 
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surface-vapor interface which could be of concern to the spacecraft 
designer. For example, these forces could overcome the small 
gravitational torques required for operation of a gravity gradient 
satellite. 

Initial work in the field of low pressure evaporation was con- 
ducted by Langmuir [ l] in 1913. He was interested in sublimation 
of incandescent light filaments and developed a semi- empirical 
expression for the rate of evaporation. Since that time many more 
experimental studies have been performed to evaluate the evaporation 
(sublimation) coefficients of a variety of materials. A compilation of 
experimental evaporation coefficients obtained through 1961 is given 
by Paul !. 2 ] . He indicated that the majority of materials which have 
been studied evaporate into a vacuum at or near the maximum rate 
given by the Knudsen-Langmuir expression for effusive flow 


m = 1 / 4 p C 


( 1 ) 


and p is the vapor density. In a vacuum environ- 
ment, the experimental evaporation coefficient is obtained by dividing 
the measured or calculated mass flux by that given in equation (1). 
Materials which evaporate at significantly lower rates than the 
maximum were generally characterized as those which exist in the 
vapor in forms different from the condensate. These conclusions 
were in agreement with an earlier study made by Knacke and 
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Stranski [3]. Paul also noted that unclean surfaces and experimental 
errors tend to reduce the evaporation coefficient so that in some cases 
the experimentally determined value may be lower than the true 
value. 

Experimental studies undertaken to determine the evaporation 
coefficient for 2-ethyl hexyl phthalate in a near vacuum environment 
led Hickman and Trcvoy [ 4 J to notice the effect of a small back 
pressure on the evaporation coefficient. They found that at a limiting 
vapor pressure of 1 p the evaporation coefficient was near unity while 
a two order of magnitude increase in the vapor pressure reduced the 
coefficient to 0.75. In a later study with water, Hickman T 5j again 
obser ved the influence of back pressure on the rate of evaporation. 

He obtained an evaporation coefficient of 0.2 5, which was consider- 
ably higher: than most values given previously. 

A number of theoretical investigations have been conducted in 
which evaporation from a surface was studied with a wide range of 
back pressures. In 1936 Crout [6] considered the one dimensional 
problem of evaporation of a monatomic vapor from a surface. He was 
primarily interested in evaluating the gas properties at the vapor- 
surface interface. To do this, Grout developed a modified Maxwellian 
distribution incorporating four constants. The values of the constants 
were determined by balancing the gas-liquid mass, momentum, and 
energy fluxes and specifying the rate of evaporation. 
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A similar analysis was performed by Schrage [ 7 ] in 1953. 
However, he assumed a different form for the dis tribution function 
near the surface. The outflow (u > 0 ) was described by a Maxwellian 
distribution function f corresponding to the emitting surface temp- 
erature. The back flow (u < 0) was assumed to be represented by 

f (1+BU) where B is related to the mass flux and U is the random 

o 

molecular velocity perpendicular to the surface. As with Croat, the 
objective of Schrage' s analysis was to evaluate gas properties at the 
interface, but the mass flux was an undetermined parameter. 

Schrage studied both monatomic and polyatomic vapors. 

More recently, Collins and Edwards [8] studied evaporation 
from a spherical surface into a vacuum or into a pure vapor under 
strong nonequilibrium conditions. The object of this investigation 
was to determine the effect of molecular back, scatter in the encom- 
passing vapor cloud on the rate of evaporation for strong nonequili- 
brium conditions. The continuum assumption was applied and both 
monatomic and diatomic vapors were considered. A Grad repre- 
sentation for the distribution function was used to permit the connec- 
tion of the surface boundary conditions and the gas dynamics in a 
consistent manner. It was found that for evaporation into a vacuum 
with infinite Reynolds number, the evaporation coefficient was 
independent of surface temperature and equal to 0.8116 for a mona- 
tomic gas a lid 0. 7778 for a diatomic gas. Evaporation into a homo- 
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geneous vapor was also studied. The evaporation coefficient was 
found to be greater than one in certain cases. In both problems the 
rate of evaporation was determined not specified as it was in Croat's 
[6] and Schrage's [ 7] case. 

Within the last two years three papers have used Lees' moment 
technique [9] to approximate evaporation and condensation phenom- 
enon. 

Patton and Springer [ 10] studied two quasi- steady problems: 
i) evaporation from a plane surface into a vapor 
ii) flow between two parallel plates at different temp- 
eratures. 

In their analysis the vapor is treated as an ideal gas composed of 
monatomic Maxwell molecules, with a Lees' representation of the 
distribution function. They employed four moments (mass, x- 
momentum, energy, and x-heat flux where x is the direction of 
motion) of the Boltzmann equation and solved a linearized form of the 
resulting equations in order to obtain an analytical representation for 
the mass flux in terms of Knudsen number for the two problems 
considered. 

In a similar manner, Sampson and Springer £ 11 ] investigated 
the evaporation of a spherical drop into a pure vapor. Again four 
moments were taken and the resulting equations linearized to obtain 
the mass flux. They also considered droplet evaporation into a eras- 
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vapor mixture. 

Both condensation of a vapor on a flat surface and evaporation 
from and condensation on a spherical drop were studied by Shankar 
[ 12]. Besides solving both problems using the quasi-steady assump- 
tion, he also obtained a solution for the nonsteady flat plate conden- 
/ 

sation problem. Shankar used the same four moments that were 
employed by Springer and his associates to determine the mass flux 
for both problems. Furthermore, he solved the liquid- vapor inter- 
face (flat surface) problem with a six moment method and found 
essentially the same expression for the mass flux as given by the 
four moment procedure. He did not consider the two plate problem. 

The advantage o.t Lees* [ y] multiple moment kinetic theory 
technique is that it affords a method of solution which is in some 
circumstances valid over the range of flow conditions from free 
molecular to continuum. The corresponding equations are, however, 
so complex that even with a one dimensional evaporation problem it 
is necessary to linearize the moment equations in order to obtain an 
analytical solution. There is, as always, a question regarding the 
range of validity of the linearization. For example, the linearized 
four moment method u3cd by Springer and his associates and Shankar 
to solve the flate plate and spherical problems does not allow for a 
near equilibrium condition to exist with mass flux. Further, the 
small parameter used to linearize the equations is the mean velocity. 
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Thus, a problem arises in the case of s umultaneous small velocity 
and high Reynolds number limits. 

In the present investigation the problem of flow between two 

pai. allel plates at different temperatures is solved numerically for 

some cases using Lees' moment method. The influence of the 
/ 

induced vapor cloud on the evaporation rate and the vapor motion 
between the two plates can be studied over a range of flow conditions. 
Six moment equations are used along with the proper boundary condi- 
tions. 

Both the non-linear and linearized problem are solved. The 
non-linear problem is a two point boundary value problem which is 
solved numerically. The solution to the two plate problem is given 
analytically for the case where the equations may be linearized as 
small deviations from equilibrium. 

Several simplifying assumptions are made. They are: the two 
plates and the vapor are composed of the same substance; the vapor 
is a monatomic gas consisting of Maxwell molecules and obeys the 
perfect gas law; the accommodation coefficient is unity, i. e. , every 
molecule that strikes the surface will be absorbed by it; and a 
Maxwellian distribution corresponding to surface temperature with 
zero mean velocity describes the molecules emitted by a surface. 

The possible effect of surface structure on the evaporation and con- 
densation rates is neglected. 
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KINETIC THEORY FORMULATION 
2.1 The Boltzman Equation 

In fluid flow fields where large gradients occur, translation 
non- equilibrium effects are observed by the presence of viscous 
stress and heat flux. For such flows the Maxwell distribution func- 
tion does not adequately represent the translational molecular 
velocity and it is necessary to obtain other expressions describing 
the physical process. Consideration of the conservation of mass, 
momentum and energy in the absence of external forces results in 
five differential equations describing a larger number of dependent 
variables. For Newtonian liquids and perfect gases at normal 
densities it is possible to empirically justify additional constitutive 
equations and an equation of state which allow the number of depen- 
dent variables in the conservation equations to be reduced to five. 

For gas flows at lower density the constitutive equations are not valid 
and it is necessary to approach the problem from a different point of 
vi ew . 

A number of methods have been developed to describe gas flows 
over a range of gas dynamic regimes. The most complete of these 
is to follow all of the particles in their collisions throughout the flow 
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field. Such attempts have met with little success because of the 
great complexity, although Monte Carlo techniques which follow 
representative molecules are being used successfully in some 
problems. In lieu of following the dynamical trajectories of separate 
particles it is mathematically feasible to represent an approximate 
variation of the particle distribution function throughout the flow. 

Such an approximate formulation is applied in this investigation and 
is described briefly below. 

The variation of the molecular distribution function is governed 
by the Boltzmann equation. This equation can be derived either by 
introducing appropriate time averages into Liouville's equation for 
an N particle system [ 9] or by writing an equation for the rate of 
change of the number of particles in a given velocity range (Yincenti 
and Kruger [ 13 ]). It has the fo rm: 


|i tc .±L + i -2- ( F.f ) 

dt j d x . Sc. j 

J J 


dA 


Lf (c.')f(cl) - f(c.)f(c.)]VdA dc„ (2) 
1 J 1 3 c J 


where f is the distribution function, V is the relative velocity between 
colliding molecules , and dA c is the generalized differential collision 
cross section. Integration is over the velocity space of the colliding 
molecules. From left to right the terms of this equation may be 
interpreted as the rate of change of the number of molecules of class 
Cj whjch results from convection, external forces, and collisions 
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with other molecules. There are two implicit limitations of the 
Boltzmann equation [ 13]. First the range of intermolecular forces 
of the gas must be small compared to molecular separation which 
must in turn be small compared to the mean free path. Such a 
limitation corresponds to the assumption of a thermally perfect gas. 
Second the distribution function must not change appreciably over a 
distance of the order of the range of interparticle forces or time 
interval of the order of the duration of a representative collision. 

For the majority of gas dynamic problems of interest these 
limitations cause no problem. The Boltzmann equation is difficult 
to solve because of the numbers of inolecules involved and because 
of the complexity introduced uy tiie non-linear collision integral. 

As a result it is usually necessary to introduce approximate methods. 
A number of techniques have been developed which exploit the possi- 
bility of linearizing the collision integral term. Unfortunately, none 
of these methods give results applicable over the range of flow 
regimes from free molecular to continuum and it is necessary to 
introduce another form of analysis. 

2.2 Maxwell's Equation of Transfer 

The difficulty inherent in an attempt to solve the Boltzmann 
equa tion directly is not the only motivation to find another kinetic 
theory formulation. Maxwell recognized that it is not the distribution 
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function itself that is of interest but certain lower moments which 
correspond to physical variables of interest. As a result he devel- 
oped an integral equation of transfer for any quantity Q which is a 
function of particle velocity. In general, such an expression may be 
derived either by considering the sources of change of Q in the physi- 
cal space or by multiplying the Boltzmann equation (2) by Q and inte- 
grating over velocity space. The resulting equation is known as 
Maxwell's equation of transfer and takes the form: 


9Q , 3 , — — , ^ 3Q ,r Al 

-t— + 3 — (c.Q) - F. = A[QJ 

3t 3x. j j 3c. 


" F J 


\\~c - ) ii j v ^ 


_ oo _ oo 


( 3 ) 


where (Q'-Q) is the change in Q(c.) resulting from a molecular 
collision and V is the relative velocity between colliding molecules. 
Integration is performed over the velocity space of both the probe 
particle of interest (unsubs cripted) and the colliding particle of a 
different class (subscripted 1). As before, dA represents the 
general expression for the differential collision cross-section. 

As with the Boltzmann equation for the distribution function the 
terms above may be interpreted from left to right as the rate of 
change of Q in a fixed volume due to particle convection, external 

i 

forces, and collisions. In this expression the distribution function 
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does not appear explicitly but rather as a weighting function, the 
form of which is discussed subsequently. In this manner the 
Boltzmann equation for the distribution function is not satisfied 
locally but rather in some average sense. Such an approach is 
analogous to the integral techniques used in solving boundary layer 
equations. 

The equation of transfer given above cannot, in general, be 
reduced any further without a knowledge of the distribution function 
and the details of the collision process. However, if Q(c i ) repre- 
sents the mass, momentum, or energy per molecule, these equations 
simplify. For such functional forms of Q conservation of mass, 
momentum, or energy during elastic impact require that the right 
hand side of equation (3) be zero. For these collisional invariants a 
system of equations is obtained which is independent of the collision 
process. These equations constitute a coupled set of five differential 
equations for the conservation of mass, momentum, and energy of a 
monatomic gas occupying a fixed volume in physical space. Depend- 
ing on the form taken for the distribution function it is generally 
necessary to consider additional moments to obtain the proper 
number of differential equations for the undetermined functions in 
the distribution function. It is in evaluating these "higher" moments 
that the details of the collision process must be specified. 
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2.3 Maxwell's Inverse Fifth Molecular Repulsion Model 

The collision integral on the right hand side of equation (3) 
may, in principle, be solved for any arbitrary distribution function 
and law of force between colliding molecules. Practically speaking 
it is desirable to choose a law of repulsion which affords the greatest 
mathematical simplicity yet retains the non-linear character of the 
collision integral and the short range interaction behavior implicit 
in the derivation of the Boltzmann equation. These considerations 
prompted Maxwell to suggest an inverse fifth power repulsion law 
which takes the form 


F 


m ^m^Kr 


(4) 


where K is a constant, r is the distance between centers of molecules, 
and nij and m^ are the molecular mass. 

The inverse fifth power repulsion law does not provide a 
particularly accurate description of intermolecular forces, however, 
it does allow one to simplify the collision integral. For this repul- 
sion law the relative speed of the colliding molecules vanishes under 
the integral and the collision integral may be written as 


where 


A[Q j = 



Jjff^JdcdCj 


J = 


oo 2rr 

i (Q'-Q)ds a da 

■ t 

o o 


(5a) 


(5b) 
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and e and a are parameters describing the collision process. Thus 
for Maxwellian molecules the collision integral may be interpreted 
as the value of J averaged over the velocity space of the two partici- 
pating classes of molecules. Furthermore, J is proportional to the 
value of Q and is independent of the velocity distribution function. 

2.4 Lees' Bimodal Velocity Distribution Function 

All tha t remains to complete the kinetic theory formulation is 
an expression for the velocity distribution function. Although the 
form of the distribution function to be used in solving the equation 
of transfer is not unique, basic requirements to be satisifed are: 

i) It must have the "two-sided" character essential to 
rarefied gas flows. 

ii) It must be capable of providing a smooth transition 
from rarefied flows to the Navier -Stokes regime, 

iii) It should lead to the simplest possible set of differen- 
tial equations and boundary conditions consistent with 
requirements i) and ii). 

Guided by the limiting solution for free molecular flow;, Lees 
[ 9 ] suggested that the distribution function take the form of a 
"two-sided" Maxwellian. At a given point the contributions of the 
two "sides" are determined by line of sight. This is illustrated 
schematically in Figure 1 for a spherical body placed in an unbounded 
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free-molecular gas with diffuse reflection at the body surface. At 
a point P, particles with velocity vectors lying within conical Region 
I are described by a Maxwellian corresponding to the velocity and 
temperature of the surface. The distribution function for the re- 
maining particles emanating from Region II is the free-stream 
Maxwellian. 

In the general case the two regions are determined by the line 
of sight principle and the distribution function takes the form 


c in Region I 


f = f = 


- r j, -, 3 / 2 [c- Ul (r,t)V 

r n, (r, t) I ~ I exp { - — — } 

v ^KTjIr.t) J 


( 6 ) 


c in Region II 


f= f 2 = 


r 1 1 3/2 f tc-u (r,t)]' 

= n 7 (r,t) — exp j- 

2 L ttRT K H J l 


2TiRT 2 (r, t) 


2RT 2 (r, t) 


} 


In these two expressions 2 (r,t), u J 2 (r,t), andTj 2 (r,t) are ten 
undetermined functions of space and time. It is then necessary to 
determine these functions by solution of ten simultaneous moment 
equations. When these functions are specified, all macroscopic 
quantities of interest can be evaluated. 



Chapter 3 


MOMENT EQUATIONS AND BOUNDARY CONDITIONS FOR 
THE FLAT PLATE EVA PORATION - CONDENSATION PROBLEM 

3. 1 Problem Definition and the Distribution Function 

The present investigation considers the two plate problem. The 
surfaces are maintained at unequal temperatures and external forces 
are ignored. Emphasis is placed on understanding the physical 
behavior of a single component vapor between the two bounding 
surfaces. It is assumed that the two surfaces are maintained at 
constant temperature. The hot wall is the evaporating surface 
(source) and the cold wall is the condensing surface (sink). It is 
further assumed that the vapor between them is monatomic, obeys 
the perfect gas law, and consists of Maxwell molecules. The process 
is quasi -s tea dy. 

The problem is illustrated schematically in Figure 2. The hot 
plate, at temperature Tj, is located at x = 0 and the cold plate, at 
temperature T , is placed at x = ■&. For this problem gradients of 
physical variables and parameters in a direction parallel to the 
plates are zero. As a result, the expression for the two-stream 
distribution function may be written 
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Figure 2 The Two Plate Problem 
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In these two expressions, n (x), u (x), and T (x) are six 

i ) Ll X y Cj 1 y Li 

undetermined functions of the independent variable x. With this form 
of the distribution function the six parameters may be evaluated by 
simultaneous solution of six independent moment equations. 


3.2 Equations of Transfer 

In obtaining these equations, it is desirable to consider the 
lowest moments of the Boltzmann equation for two reasons: first, 
they allow the greatest degree of mathematical simplicity; and 
second, the "lower" equations generally involve moments of the 
distribution function which may be interpreted as physical variables 
of interest. With this in mind, the six lowest independent moment 
equations will be developed; this will involve choosing successive 
fo rms for Q(c^) that represent the lowest powers and combinations 
of molecular velocities. 

As discussed previously, choosing Q(c.) to be the molecular 
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mass, momentum, and energy results in moment equations which are 
conservation equations for the respective quantities. For these 
values of Q(c.) the collisional term is zero and in the absence of 
external forces the generalized moment equation (3) becomes 
successively 


^(p) + a ~-(p c ) = 0 

J 


( 8 ) 


¥ (pc i ,+ sT (p c j c i> = 0 

J 


(9) 


a 2 S 2 

3F l ' c l , k l,c j t| = o 

J 


( 10 ) 


2 --v 2 

U T r* ~ ^ i ~ - J-*., -V-.rl .... .-1.1 -» , r-C 

— «— v- . i j_ oi. tile ottuuv - Otalh OuC Uiiiiciio j.uilat x c/ »_> a. v... jl ». j. ui 

i= 1 1 

interest, time and the y- and z-derivatives are zero. Thus, the first 
three independent moment equations are: 


Q 1 = m: 


— (pc ) - 0 

dx x 


Q„ = me : 
2 x 


£ (P C x> = 0 


= nnc /2: 


d , 2 V „ 

d^ (P V > = ° * 


( 11 ) 


( 12 ) 


(13) 


Expressing these equations in terms of physical variables is 
facilitated by equating the molecular velocity to the sum of the mean 
and thermal velocities 
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c=c+C=u+C. (14) 

x x x x 

In terms of the thermal velocity the following moments of the distri- 
bution function can be identified: 


a.. = -p C.C. = t . . - p 6 . . 
U i J ij ij 


( 15 ) 


3 

2 


kT 



(16) 





(17) 


In addition, the following equations are appropriate for the distribu- 
tion function and gas model assumed 

n = J f dc (18) 

p = nkT = p RT (19) 

h = e + — = ^ RT . (20) 

p 2 


By using equations (13) to (20), the first three moment equations may 
be rewritten as 


dx 


[ p u 3 = o 


d r 2 no 

-T- L PU - CT 3=0 
dx xx 


h[p U (| RT + = 


( 21 ) 

( 22 ) 

(23) 
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Integration of these equations yields 


p u = mB. 


2 

pu - a 


= mB. 
xx 2 


- 4 . 

S 3 U X 

>u ( — RT + — j + q - o u - mB 
V 2 2 y x xx 


(24) 
(2 5) 
(26) 


where B , B > and B are integration constants. 

X. Cm ' J 

The next higher moment equations are obtained by setting 

Q Ac.) = mc.c, . 

4 l j k 


The resulting equation of transfer is interpreted as an expression for 

the flux of momentum, me, , in the x. direction. For this value of 

k J 

Q the collision term on the right hand side of equation (3) is non-zero 
and is evaluated by Lees [ 9 ] for Maxwell molecules. In this 
problem, equation (3) simplifies to 


d r 

— — |_ m c c .c 
dx x j 


k 


] 


— T .. 

H jk 


(2 7) 


where j and k may independently assume values from 1 to 3. Three 
non- trivial equations are obtained for j = k, although due to symmetry 
in the y and z directions the two independent equations resulting are 



2 

me : 


x 


_d_ 

dx 


C me 3 ] 
x 


P 

U 


T 

XX 


(28) 
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Q = me : 
4b y 


d r 2 1 _ P 

-t— L m c c J - — t 
dx x y (J yy 


(2 9) 


Adding twice equation (29) to equation (28) gives 


d 

dx 


r , 2 2 2 .i 


me (c + c + c ) 

= - ( 

xx y z _ 

M V 


l Py' = 0 


which is identical to the energy equation (13). Thus of equations (13), 
(28), and (29), only two are independent; the choice of the two to be 
used in the solution for the arbitrary functions will be made on the 
basis of simplicity. 

The fifth equation of transfer is obtained by setting Q (c ) = 

4 i 

2 

o. *~ (<- / ^ ) uau uio. j( is** .ubipn.Cu uC the * 1 ux o. energy m the j- 

J 

direction. Lees [ 9 ] has evaluated the corresponding collision 
term for Maxwellian molecules and in this case the resulting moment 
equation becomes 


Q = m c (c /2) : 
5 x 


_d 

dx 


r 2 c ■ 

1 P 

2 -| 

-me — 

= - 

- — q + T u 

L x 2 

J M 

. 3 M x xx J 


(3 0) 


The sixth independent equation of transfer is obtained by 
3 

choosing Q^(cg) ~ mc j • The meaning of this equation is less physical, 
as is typical of higher moments. It represents the flux of me 2 in the 

j 

j-direction. The collision integral is evaluated in Appendix A and 
the moment equation reduces to 



24 


^ 3 

Q, = me : 

6 x 

In summary, the six independent differential equations of trans 
fer to be used are 


_d_ 

dx 


r 4 i 

= - f 

q + 3 ut - ~~ m 

cH 

L x J 

M 1 

. x xx 2 

X J 


(31) 



dr 3 1 p 

- — me | — i 

dx L x J (i xx 


A 


u 


two of 
three 


yco, 


iL r 

dx L 


2 

m c c 
x y 




— T = 0 
xx 


J 


(2 9) 



(30) 


(31) 


3,3 Expression of Physical Variables and Moment Equations in 
Terms of Parameters of the Distribution Function 
In order to solve the above equations for the parameters ri^, n^ 
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U 2 ’ j * ^2 functions of x, it is necessary to express the 
equations in terms of these variables. In principle this involves 
nothing more than evaluating the prescribed moments of equation (7) 
and substituting them into the equations of transfer. In performing 
these operations and in reducing the moment equations to their 
simplest form, a significant amount of algebraic manipulation is 
required which adds nothing to the understanding of the problem. 
Therefore, the contributing physical variables and higher moments 
of the distribution function are evaluated and the final forms of the 

equations of transfer are presented without details of the intermediate 
algebraic steps. 


7he 


appiopxiate uj.oi.u6i.ils of Ine Liuioual distribution function 


(7) were obtained by using integrals summarized in Appendix B. 
Because of the recurrence of certain functional forms in all moments 
it is convenient to define 


E l" 

1 + erf 

r u i ( x ) -| 

L ^2 R T 1 (x ) J 

(32a) 

E 2" 

1 + erf 

- _ u 2 ( x ) n 

_ /2RT (x) J 

(32b) 

X ! S 

exp -u^(x)/2RT (x) 

(32c) 

X 2 S 

exp -u 

|(x)/2RT 2 (x) 

(32d) 
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and to replace the variables T^(x) and T^(x) with 

= 2RTj(x) (32e) 

P 2 = 2RT 2 (x) . (32f) 


With these definitions, the physical variables and higher moments 
may be expressed as 


p -* 1 1 

n - J f dc = — n , E , + — n„E, 


2 1“1 2 22 


nu 


= J fc , d « ■ K(“i e i + /? x i) + .K( w/t x z> b i 


a = -m 

XX 


fC 2 dc 

x 


-m 


f(c - u) dc 
x 


-J mn z{[-i + 4] s z- a zl^ x z} 


+ mnu 


p Z 

o = a = - m f C dc 

yy zz J y 


= -imn 1 B 1 E 1 -imn 2 6 2 E 2 


P = " 


i ( o u /s, = 
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imn 2 [( u ^|p 2 )E 2 - u 2 /Jx 2 ]-I 


mnu 


. m p _ 2, 

q - - C C f dc 
x 2 J x 


X 


- 4"""! { - [| (u-u 1 )Bj+(u-u 1 ) 3 ]e i + [2S 1 +Uj + 3(u 2 -uu 1 )] Jjj. 

~ 7 mn 2 {[ ! W < u -“ 2 > 3 ] E 2 + [ 2 B 2 + -.2 + 3 („ 2 -uu 2 )] yfxj 
J £c * d ‘ = I“i{( u i + I u i b i>i + (Pi +u i)/t" x i} 

+ I n 2 {( u 2 + 1 u 2 s 2 ) E r ( V U D /r x 2 } 


["i E . + / 5 X l] + r' 2 B 2 [“ 2 E 2 -y¥ X 2 ] 


, 2 ,- 1 

f c c dc = — n p 
x y 4 11 


f c 


2 c 

x 2 


1 r r 5 2 2 4~i r 3 7 "l/^l i 

dC = 4 n l{l4 S l + 4 Vl + u l] E l + [“, + 2 u l S lJ •/— X 1 } 


1 rr 5 2, , 2 4*1 r 3 7 

4 n 2 i|_ 4 B 2 +4 "2 B 2 + u 2 J E 2 - [ U 2 + 2"2 B ; 


•A _ 
TT 2 


} 


4 -* 
f c dc 
x 2 


|"l{[ l e i +3e i u l + ^] E l + [ u i + |u 1 B 1 ]/4 Xj} 


+ 


? "2 {[ I S 2 + 3S 2 U 2 + 4 ] V [ u 2 + I " 2 B 2 ] /i” X 2 


} 
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J f dc = J f(c^- u)dc 

= i n i{ - [l e i +, "- u i> 2 ] (u “ u i )E i +(B i +u i +3u2 ' 3uu 1 ) /ir x i} 

+ \ n 2 {'[ I B 2 +(u - u 2 |2 ] (u ‘ U 2 )E 2 ' (B 2 + “2 + 3u2 - 3u “ 2 > / 


“ X 
TT 2 


} 


Therefore, the three conservation equations become 


d 

dx 


>i( u i E i + 





i[ 2B i] = 0 < 33 > 


d 

dx 





+ u. 






} 




U 2 S 2 


+ u 2> : 




[ 4B 3>' 


(3 5) 


Similarly, the remaining moment equations may be simplified by 
combining terms and by utilizing the integration constants B^, 

and B^ to yield 



ab[ n i (u M u i B i )E i +n i (e i +u i ) 



T= I'D 
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+ >’ 2 C'2 + | " 2 B 2>V n 2<V U 2> j T x 2 ] 

-|.£[2„u 2 -2B 2+ =0 


+ I M L 2 " uZ " 2B 2 + ? (n l E l V n 2 E 2 6 2 ) ] = 0 


E[ n i(! s i + 4 u i B i +u D E i +n i( u i + ? u i e i)/^ x , 

+ n 2 ( ! B 2 + 4 u 2 S 2 + U 2 ) E 2 - "2 ( U 2 + \ “ 2 P 2 ) /v X 2 ] 

+ f m[ 2B 3-“< B 2 + Vi E 1 + “2 B 2 E 2>] = ° 


d 

dx 




+ n z(f B 2 +3P 2 U 2 +u 2) E 2- n z( U 2 + f U 2 B 2 )/t X 2 ] 


(36) 


(3 7) 


(38) 


(39) 


{-u[B 2+ i(n 1 B 1 Ei + „ 2 B 2 E 2)] + 2 B 3 -|n 1 B 1 (u 1 E 1+ y^x 1 ) 
" 4 n 2 S 2 C U 2 E 2‘ fV X 2^ } = °‘ 
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3.4 Non-Dimensional Equations and Boundary Conditions 

Non-dim ensionalization of the physical variables and moment 
equations is accomplished by introducing the following definitions 


X = 

d x 

n, ^ 

= n n, „ 

1,2 

I 1,2 

u_ __ 

= U u, , 

(M 

r— 1 

1,2 


2— 


= U P lf2 


where the characteris tic velocity U = /RTj and variables subscripted 
I are evaluated at the surface of the hot plate. Note the cold plate 
could just as easily be used as a reference. The viscosity law for 
a gas composed of Maxwell molecules (Appendix A, equation (A14)) 
is 


1 i = 



With this result, the coefficient 


P 


d 

U 


becomes 


p d _ n Re 

m” u /ym 


n Re 


where the Reynolds number is defined by 


Re = 


Pj Ud 


and /"v M = 1 because M = 


U 


/ Y RT t 


( 40 ) 
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For spherical molecules the relation between viscosity and the 
mean free path X is 

1 - 

|J = j p C 1 - 

Hence, the Reynolds number in this problem can be related to 
Knuds en number (Kn) by 



(41) 


where Kn = X/d. 

In terms of the non-dimensional quantities, the integrated 
conservation equations and physical variables of interest become 


continuity 



u iV 



} 


+ 



{ j 2 E 2 



(42) 


x-m omen turn 
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+ T'[C! U 2 3 2 +U 2 3 ) E 2"O 0 2 +U 2 2 )'/-| X 2 ] 


n = I n i E l + \ n 2 E 2 P/T 


pressure 


p = i"i [ C u f + r B 0 E i + “i 7 v x i] 


+ | n 2[( U 2 + l P 2) E 2- U 2Vv X 2 ] • T" u 


1 --2 


normal shear stress 


3 n l[ U 1 E 1 + J7 X J ' 3 n 2[ U Z E Z~ U Z r f^ X z] + 3 n U 


2 2 


heat flux 


q x b ! n i{-[l (u - u i ,B i + (u - u i ) 3 ] E i + [ 2 V u i + 3 <u 2 - uu i ) ]/^' x i}. 


"? n 2{[l (u_u 2 1B 2 + (u ‘ U 2 )3 ] E 2 + [2S 2 +u 2 + 3(u 2 -u u,,) 


~ X 2 


The evaporation coefficient is given by 
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2 /tt B 


1 1 1 / UX ' A 
4 P IV TT ‘ 4 P II V TT 


The remaining four moment equations are 


x-m omen turn flux 


n l[ U M u l B l] E l + ‘' 1 (e i +u l 2 )/¥ X l + " 2 [“2 + l“ 2 B2]E : 


n 2 (B 2 + u 2 2 , 


P 2 2 

"tt~ X 2 } - jRe n {2n u 2 - 2B.,+ - (n^Ej + n 


y-m oin en tnm fl nv 


® [ "i B i( u i E i + V V x i) + " 2 e 2 ( W •/ V x 2 )] 


2 -r 2 - 1 _ 

3 Ren [ 2 nu _ 2b 2 + I (n i p i E i +n 


2 f3 2 E 2 ) J ° 


x- energy flux 


dS [ n l( ! 9 1 Z+ 4 u 1 2 5i + “ 1 4 ) E 1 + '’i(-> 1 3 +|u 1 6 1 ) 7 -^X 1 


n 2 (! B 2 + 4 “ 2 V u z)V^( u Z + r“ 


n 2 ( u 2 3 + l u 2 6 2 )v V x 2 


| Ren [ 2 B 3 - u{B 2 +n 1 B 1 E 1 +„ 2 B 2 E 2 )] = 0 
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x-flux of c 


d r-/'3— 2 „ 2 -4\ - r— 3 5 

df ["iC? 6 1 + 3B 1 U 1 + U 1 ) E 1 + n i( U l + 2 U 1 B 1 


- — 3 . 5- 



-X 

TT 1 


+ ‘'2(l B 2 + 3e 2 U 2 + U 2) E 2 - n 2 ( U 2 + | U 2 P 2)V TT X 2 ] 


3.5-- 


+ 2 Re n { - u[ B 2 + ? (n^E,* n^E.,)] + 2B ; 


-| n l S l( U l E l + VV X l)-| n 2 8 2( U 2 E 2-yV X 2)} =0 ‘ ,53) ■ 

The boundary conditions for this problem are established direct- 
ly from the assumptions related to the nature of the flow at the two 
walls. In the absence of any conclusive evidence to the contrary it 
is assumed that ail incident molecules are absorbed and thermally 
accommodated at the wall surface. Furthermore, molecules are 
emitted from the surfaces with zero mean velocity and with a local 
Maxwellian distribution corresponding to the plate temperature and 
the number density of the saturated vapor at that temperature. 

Mathematically these assumptions are equivalent to the 
following boundary conditions 

at x = 0 


n l =n i 


n i = 1 


v° 


U], = 0 


(54a) 
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at x = 1 


and T i * T n- 


T = T 
1 I 


B i = 2 


n 2 n !I 


n 2 = W 


u 2 = o 


u 2 = 0 


T = T 
2 II 


e 2 = 2t u / Ti 


(54b) 



Chapter 4 


THE LINEARIZED TWO PLATE PROBLEM 


In order to find an analytical solution to the non-linear 
equations (42) to (44) and (50) to (53), subject to the boundary condi- 
tions (54a) and (54b), a small deviation from equilibrium is consi- 
dered. This approach, and the resulting equations are identical to 
those given by Shankar [ 12 ] but different boundary conditions will 
be applied. The first order perturbation solution is found by using 




u 

1 


u_ 


= 1 + N e + . . . 
= 1 + N e + . . . 
- i j. 7 ** - -j- 

• L r • 

= 1 + t e + . . . 

= eUjf ... 

= e U + . . . 



where e is a small parameter and N , t , and U are of 0(1). 

A | ^ 1 | u I| L 

Definitions (32a) through (32f) become 


Pj = 2(1 + t e + ...) 


V 

2(1 + t € + . . . ) 


E 1 = 

1 + V 

^ne + ... 

TT 1 

(56) 

E 2 = 

‘-V 

1 - U, e + . . . 
TT 2 - 
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Xj - 1 + 0(e ) + . . . 
X 2 = 1 + 0(e 2 ) + . .. 


where 3 - 2T. The boundary conditions are 


at x = 0 


Nj (0) = 


t r (0) = 


<n H' n l ) 

AN 

"II 

n H 

< T H- Tj, 

AT 

T II 

~ T 

II 


AN 


= - AT 


U^O) 


= 0 


at x = 1 


iyi) 

•t 2 (i) 

u 2 d) 


= 0 
= 0 
= 0 


(57a) 


(57b) 


where AN and AT are of 0(e). The cold wall conditions are used as 
the references values in solving the linearized problem. This is 
consistent with Springer and Patton [ 10 ] . 

By placing equations (55) and (56) into the moment equations 
(42) to (44) and (50) through (53), we obtain 


2<N r n 2 > + T r V 2 7f <v u 2 > = B f 


(58) 


Nj+ N 2 + Tj+ T 2 + 2 


<V U 2 ) = 2B- 


(59) 
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N 


r N 2+ 3/2(T r T 2 ) + 5/4^- <V U.,) = B> 


n — 


(60) 


h (N r N z + 3/2(T r t 2 ) + 3/2 Vr < V u 2 » 


■ -VT 


< u r u 2 > 


3 Kn 


(61) 


^(N r N 2+ 3/2(T r T 2 ) + 7f (U 1+ U 2 )) 


/— (U i " U ) 

= i/s/- -i L. 

' V 2 Kn 


(62) 


4 < 5 < N i + V + 10(T 1 + V + 12 /F {u r u 2 » 


(63) 


i r 

3 Kn 


i>V V ■ 7,T r t 2>] 


i[(3(N 1+ N 2 ) + 6(T 1+ T 2 ) + 8 /I (Uj- U 2 ))] 

"zk ( 2 ' N r N 2>- < T r T 2>) 


(64) 


where 


B 1 = 4 Vf B 1 


B, 


B 3 


1 +B ' 


2 B 3 



Re 


■fi 


IL * 

2 Kn 


(41) 
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and Kn = \/d. Define 


N + - Nj+ N 2 N_ = Nj- N 2 

fc + = fc l +t 2 L = T 1~ T Z (6 5) 

V + = v x +7 2 v_ = v r v 2 

where lh - /2tt v.. As previously indicated, only two of the three 
equations (60) to (62) are independent. The six equations that 
Shankar [12] used to study condensation at a liquid-vapor interface 

are utilized in this work. With the aid of the above definitions, we 
have 


2N_ + t + 2tt = B',' 

“ ' A 

(66) 

N + + t + + 4 V _ = b" 2 

(67) 

it - 2N = 2B" 
3 

(68) 

^(2N_ +3? + 2„v + ) = 2/3 UL- 7 

(69) 

d b" 

5= (5N + 107 + 24 V ) = .2/3 3 

“ Kn 

(70) 

~ (3N + + 6t" + + I6v_) = A (3 T_ - B” ) 

(71) 


where equation (68) is found by combining equations (58) and (60). It 
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should be noted that equation (68) is exactly = constant where 
q (heat flux) = q^e + . 

The solution to the above equations is 


t = 

^ -Ax , _ +Ax , 3 

De + Ee + — 

b 

(72a) 

V 

„ , n n~ , r n~ +a * 

“l ' 15Kn x + 3D V 5n e ' 3E V 5? e 

(72b) 

V = 


(72c) 

V 

D 11 T3 * * _ _ 

1 3 4D -Ax 4E +A x 

2tt 5rr n tt 

(72d) 

N = 

7/2 D e ‘ Ax + 7/2 E e* Ax - 3/10 B” 

(72e) 

N + = 

2B x r-r- r-y— ^ 

BL - »,t 3 +3D(2 ,/, 5 - J] >' 

2 1 15Kn V V 2iT V 5tt / 

(72f) 





where A= J —— and D, E, and or, are constants to be deter- 

2Kn V 2 1 


mined. By using the definitions given in (65), we find that equations 


(72a) to (72f) can be written as 
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2B” x 

2N 1 = B ;. Bi .3/10BJ + T5E + 


0' 2 + =:> 


■Ax 


/mr 


(73a) 


( 7/2 - -H-.') 
V /TOtt ^ 


E e 


2N„ = 


- 2B 3 X x 

B 2- °'l + 3/10 B 3 + TTS + ( 


24 


IOtt 


7/2^)d 


-Ax 


(73b) 


-0« + :=) 

X / 1 A ' 


E e 


+ A : 


/ IOtt 


2 b! 1 X 


2t r 


O' 


1 1 5Kn 


3 “ + fl + -$-> e' Ax +(l-J- >e Ax + ^ 
v -^To" y v /Tott ^ 5 


(73c) 


2t _ = 


2B 3 X B 3 + 
®1 15 Kn " 5 


B" B," 


(—-l) 

v /To 7 T J 


-Ax 


2 V = ~ + — - (- + 3/2 
1 2 rr 5 tt V tt 


■fl') 


De *■“ - (l + — ^~')Ee Ax 

V /TOW ' 


(73d) 


D e 


_ B 1 . B 3 fA 3 [Tn ■ 
2V 2 2tt + 5n " Crr" 2 V 2ir/ De 


(• 



(73 e) 


(73f) 


Application of the boundary equations (57a) and (57b) to equations 
(73a) to (73f) yields 


-it , s 7 24 \ z' 7 24 "N 

' (l + ^-) D+ ( 2 -^=r> < 74a > 


-2AN = - ffj- 3/10 B 3 + + 


7 24 


/TOtt 


/ IOtt 
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2B 3 
15 Kn 

(74b) 

(74c) 


(74d) 


(74e) 

(74f) 


These six equations can be solved for the unknown constants B^', B^» 
b”, D, E; and a. 

Let D = D'e^ and E = E'e Then: 


AN 


(- 


b: 


24 

/TOtt 


A 

sinh — - 4 cosh 


96 


V 


/TOrT 


sinh — + 14 cosh 


£> 


-132 
5 /Tott 


sinh — - 4 cosh j 


8 sinh A/ 2 


/TOtt 


Kn 


— cos h A/2 
5 Kn 


(7 5a) 


D' = E' 



2 

15 


■k> N + 


(! 


15 


1 

Kn 


) 


AT 


-132 . . A A 8 

sinh — - 4 cosh — - 

5 /IOtt v/TOrr 


sinh A/ 2 
Kn 


6 cosh A/ 2 
5 Kn 


(7 5b) 
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a . 


■a 


— 1 ^ b" 

15 


) B" + 2D' cosh ^ sinh A 


/IOtt 


B 


b' 1 = - + 9D' cosh j + sinh £ 

10 2 /ToW 2 


(75c) 


(7 5d) 


, m _ 


2 - 


A 


- — B' 1 + 16D' cosh — + 6tt D' sinh — 


5 3 



A 


C TT 


(7 5e) 


Now, it is possible to examine the two limiting extremes Kn -*0 
and Kn -* 03 . For Kn -*0, we have 


B = -1. 66 A P , AP = (AN + AT) 

V -f-i' 83 ' lp ' 


B 2 = 


-AP 


B' - 0 


(76a) 

(76b) 

(76c) 


and for Kn -*“> , we find that 


B^ = -2AN - AT 

®1 = “i/f " AT) (77a) 


B^ = -A P (77b) 

B ^ = AN - | AT . 


(77c) 
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When Kn -* 00 (free molecular limit), the mass flux to 0(e) agrees 
exactly with the Knudsen-Langmuir expression. In the continuum 
limit Kn -*0 the mass flux to 0(e) becomes independent of the temp- 
erature difference and depends only on the pressure difference. 

This result agrees in form with Fuchs [ 14 ] and Shankar [ 12 ] . 

Finally, the behavior of the density and velocity in the flow 
field between the two plates is established. The number density is 
defined by 

n = l/2(n 1 E 1 + n^) . (45) 

In terms of perturbed quantities, (4 5) becomes 

P~ = 1 + (1/2 N + + V )c + . . . (78) 

where p = mn. From equations (72c) and 72f), one finds 


P = 1 + 


Bl a. 

r Jl _l 

B. x 

+ 3 + / 

- 12 

- /3m 

L 2 2 

1 5 Kn + V 

v /TotT 

2 V 2tt J 


D'e 


r 10rr 




(79) 


i r-. , A/2 -A/2 

where D = D'e , E = E e 


AN 


D' = E' = 


(MM M-(MM) 


5/T0tt 


132 . A . . A 

sinh — - 4 cosh — 

L* 4 


8 sinh A/ 2 Jb cosh A/2 

Kn “li Kn 

(7 5b) 


/lOn 
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and and o' are given by (7 5c) and (7 5d) respectively. 
(79) at the cold wall can be reduced to 

P_= 1 + (-|o § 3 + i D ’ COSh I + 3 /^ D ' =inh|)c + 

At x = i (cold wall), we find 
Kn - 0 


P -1 - 


-Lfl 

15 V 2 


+ AP + 

/TOtt ^ 


(~t) 


/To 


TI 


Kn 


P = 1 - - AN + 


Similarly, for x - 0 (hot wall), there results 




20 15 Kn 


- 3 


X + (f 


. A , 48 . . i 

cosh — r — ■ — sinh - 

^ /TT~ & 


/I- l) D '] 


/Uh7 


e + 


Equation 


(80) 


(81a) 


(81b) 


and the limits are 



46 


or 


Kn - 0 


P = 1 - 


AN^ll + — 


IOtt 


>Af(7 + -^- 

/To^ 


:) 


15 






IOtt 


, A N (2 7. 07) + AT (12. 36) 

1 " 39. 42 


(82a) 


Kn - “ 


p = 1 - j AN + 


(82b) 


The mass velocity is 


u 


__ - fc 7 ' u l - bz _u 2 /B 2 

n l U l E I +n 2 U 2 E 2 +n l V~ e - n 2 -J— e 

"lV" 2 E 2 


(83) 


and in terms of the perturbed quantities u becomes 


- 1 [2 - , 1 - J - , , 

U = 2 VT (N - + 2 l - + " v + )e + 


(84) 


where N , t , and v + are given by (72e), (72a), and (72d) respective- 
ly. At any point in the flow field between the two plates and at x = 0 


or x= 1, (84) reduces to 
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- l [T B i - 

u= 2 Vn • — e+ '“ =B l e +, ‘*- 


(8 5) 


The limits for B 1 as Kn -*0 and Kn - 00 are given by (76a) and (77a). 

The density correction to 0(e) at the hot wall, equation (82a), 
depends. on the relative size of AN versus AT. Both AN and AT are 

negative by definition. If the density correction to 0(e) is greater 

C at" 

< 2.19); but if 

the correction to 0(e) is less than zero, the flow is temperature 

, . A at 

dominated ( “”=r > 2.19). 

VAN 

From equation (8 5), it is apparent that the small parameter e 
ife related to the mean velocity u which results from a small deviation 
from equilibrium. The fact that u is constant to 0(e) is not surpris- 


ing since 


p - = i t.r* 1 * + ,%-«*>+.. . 

u = C u< 1 > + e 2 u< 2 > + ... 


and pu - const. Therefore, we have 


; u ( 1 ) + e 2 (p< 1 >u' 1 l + u« 2 )) + --- 


const. 


u' 1 ' . c 


1 


u< 2 > = c 2 - 


which implies 
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but p ^ t constant (as previously demonstrated, Kn -* 0) so 
“( 2 ) , 

u T constant. 

The characteristic Re, Equation (40), is based on a velocity 
U = / RT^ which is of 0(a) where a is the speed of sound. Re can be 
written as 



where v = ratio of specific heats, M, - Mach number, Re is 

11 pu 

Reynolds number based on the mass flux pu at the cold wall and the 
mean velocity u is small. For extremely small kinematic viscosi- 
ties v, large d, or a combination of both, Re could be large when 
u is small. 

When u is not small, the full non-linear equations must be 
s olved. 



Chapter 5 

NUMERICAL SOLUTION TECHNIQUE 

5. 1 Solution of Separated Boundary Value Problem as Initial 
Value Problem 

The six moment equations to be solved for the unknown 
parameters n^ ^ ^ u j 2 form a system of first order ordinary 

non-linear differential equations. Because all dependent variables 
are present in each of the equations they are completely coupled and 
must be solved simultaneously. Although three of the equations are 
integrable, no method was devised to simplify the solution by using 
the mixed algebraic and differential equations. Attempts to reduce 
the number of dependent variables by the introduction of groupings 
were unsuccessful and no feasible analytic integration technique was 
developed. 

Numerical solution of such a system of equations can be 
accomplished by a number of methods if the values of all the varia- 
bles are prescribed at one of the boundaries, i.e. , if it is an initial 
value problem. With the present problem three of the boundary 
conditions are given at each surface and as such represent a 
separated or two point boundary value problem. Solution of this 
type of problem generally involves either quasi-linearization as used 
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in the last chapter or reduction to an initial value problem. Because 
the free molecular limit allows an obvious choice of initial values 
which systematically vary at higher Reynolds number, the latter 
approach is used in order to solve the problem for strong non- 
equilibrium conditions. 

The numerical solution begins at the first plate (x=0). Here, 
values of n^, u^ and 3^ are prescribed and the values of the other 
three variables n^» u^ and 3^ are assumed. Next, the derivatives of 
the six functions are evaluated at the wall. With these derivatives 
known, the six equations can then be simultaneously integrated to x+s 
by using a fourth order Runge-Kutta scheme. Here s is the normal- 
ized step size. From this point the process of evaluating the deriva- 
tives and of integration is continued to x = 1 . Atx = l the integrated 
values of n^> u^ and are compared to the actual boundary conditions 
prescribed for the problem. If the integrated values are not suffi- 
ciently close to the boundary conditions new values of n , u , and B 

L* L* Ct 

at x = 0 are calculated and the equations are integrated again. The 
iteration for the correct values of n^, u^, and B >2 at x =0 continues in 
this manner until the boundary conditions specified for n^, u^ and 
3^ at x = 1 are satisfied. Since the six parameters are evaluated at 
each step of the integration by using the definitions in the previous 
chapter, the values of the physical variables of interest are deter- 
mined throughout the flow field. The specific details of the 
integration method are briefly outlined in the next two sections. 
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The complete computer program is given in Appendix C. 


5.2 Evaluation of Local Derivatives 

To simplify the mathematics of the problem it is convenient to 
define the function Y which represents the six dependent variables so 
that 


Yj = n^x) 
Y 2 " 

Y 3 " 8 >> 
Y 4 “ 

y 5 = u 2 ( x ) 
Y 6 H 8 2 {x) * 


(87) 


Furthermore, it is desirable to leave the six independent moment 
equations in the differential form 


F i/ Y j )W j + G i (Y j’ Re ^ " ° ’ i = j = 1, 6 (Appendix C) (88) 


where the derivative W is 

j 


dY. 


dx 


1 = W. 


j = 1,6 


and F and G. are algebraic functions only of the six variables Y 

j 

and the Reynolds number. 

The subroutine DRVTV (Y.W.Re) is used to evaluate the local 
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derivative W for input values of Y. and Re (Appendix C). To do this 
DRVTV calls another subroutine DSIMQ which solves the six equa- 
tions (88) simultaneously. 

5.3 Fourth Order Runge-Kutta Integration 

A number of numerical integration techniques were considered; 
however, because of the simplicity in the present application a fourth 
order Runge-Kutta integration scheme is employed. For a system of 
two first order ordinary differential equations 

y' = g^x.y.z) , z' =g 2 (x,y,z) 

the fourth order Runge-Kutta integration for step size s is given by: 
y n+l = y n + |-(k 1 +2k 2 +2k 3 +k 4 ) +0(S 5 ) , 

Z n+1 = Z n + i ( V 2 V 2 V V +0(s5)l 


where, 


k , - s gi ( x »y » z )» 

i i n n n 


f, = sg (x ,y , z ), 
i c n n n 


k 2 = sg i ( Vi s - y „ + l k i- Vl'l 1 


1 


*2= Sg 2 ( V 2 S - y n + 2 k l’ z „ + 2 V' 
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k 3 = Sg l ( V I 3 ’ y „ + ? k 2' V?V’ 


h = sg 2 <x „ + i s - ' , „ + r k 2 - VtV' 


k 4 = sg l ( V s - W V V- 


- s g 2 (x n + S * y ~ + k ’’ z - + 


n 3’ n 3' 


This system of equations is now generalized to allow integration of 
six simultaneous equations. 

Equation (88) can be rearranged to 

Y j = -< F " 1) ij G i = g i«*. Y j } i = j = l,6 (89) 

where the value of g. is obtained simply by calling DRVTV (Y,g,Re). 

Y. at x+ s is found from 
x 

Y^x+s) = Y.(x) + 1 (k. 1+ 2(k u + k. 3 ) + k. 4 ) + 0(s 5 ) (90) 

where 


k il= sg.fx.Y.) 

(91a) 

k i2 = S Si<* + ^ s - 

(91b) 

k i3 sg i<* + I s ’ Y j + | k j2> 

(91c) 

k i4 = sg i^ x + s * Y j + k j3^ i = j = 1, 6. 

(91 d ) 


To obtain explicit values of the k's in the program, use is made of 
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the flexibility of the subroutine DRVTV. At point x values Y. and 

J 

W (Y ) are known, hence 
. i J 


k 


il 


sW (Y 


.) • 
J 


By defining 


z. 

1 


Y. + 

l 


2 



the subroutine DRVTV (z,W, Re) is called to obtain 


w i' Y j + ? k ji> ■ «i<V ?v 

hence 

k._ = sW.(Y.+ i k ) . 
i2 i j 2 Jl 


Similarly, the remaining terms k and are obtained by success 
ively setting z. equal toY.+ 4 k. _ and Y . + k, and calling for 
W.(Y.+ k^) and W.(Y.+ k.^). Once these are known Y^(x+s) is 

obtained from equation (90). The step size is chosen to be 1/Re or 
0.Q1, which ever is smaller. 



Chapter 6 


RESULTS OF NUMERICAL SOLUTION TO THE 
NON-LINEAR TWO PLATE PROBLEM 

In this chapter the two plate problem is solved in certain cases 
for strong non- equilibrium conditions by using the numerical integra- 
tion procedure discussed in the last section. Three cases are 
considered: 


Case I, 

T n (d) = 

1/2 T^O), n n (d) = 

1/2 n (0) 

(92a) 

Case II, 

T u (d > = 

1/10 T^O), n n (d) = 

1/2 n (0) 

(92b) 

Case III, 

V d > ■ 

1/2 TjfO), n n (d) = 

1/10 n^O) 

(92c) 


For all cases , u^O) = u^d) = 0. Case I is dealt with in depth 
whereas the other two cases are discussed only to illustrate certain 
differences that occur at low Reynolds number. For large Re, all 

cases considered behave in the same manner. The reference values 
are those at the hot wall. 

The direction of integration is opposite to that of the vapor 
flow. This is analogous to the situation that exists in numerically 
solving the one dimensional viscous shock equations, VonMises [15], 
In the shock problem the direction of flow is supersonic to subsonic. 
Because a saddle point exists on the subsonic side of the shock and a 
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nodal point on the supersonic side, the direction of integration is 
always taken to be subsonic to supersonic. An attempt was made to 
determine if such singularities existed at the hot and cold walls but 
due to the complex nature of the equations no definite results were 
obtained. Despite this, the numerical procedure employed to solve 
cases I, II and III simply will not march forward from the hot wall 
without eventually blowing up at some point in the flow field for 
Re > 2. Therefore, the cold plate is positioned at x = 0. 

The six moment equations used in the numerical analysis are 
equations (42), (43), (50), (51), (52), and (53). These equations 
along with the expressions for density (45), pressure (46), t (47), 
heat flux (48), and the evaporation coefficient (49) give a complete 
picture of the flow field subject to the boundary conditions pres- 
cribed (54). 

As pointed out before, each case is solved as an initial value 
problem. A one dimensional array, ALF, contains six elements 
corresponding to n^O), u^O), 6^0), n^l), u 2 (l), and 8^1). 

These are the boundary conditions specified by (54) and either (92a), 
(92b), or (92c) depending on the case considered. Note that the 
subscript (1) and (2) are reversed when the cold wall becomes x = 0. 
YO(4), YO(5), and YO(6) are the initial values guessed at x = 0 for 
n 2 (0), u 2 (0), and e 2 (0). WithALF(l), ALF(2), ALF (3 ), YO(4), 
YO(5), and YO(6) given and a Re specified, the numerical integration 
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proceeds forward from x = 0 to x = 1 in steps of 0. 01 by using a 
fourth order Runge-Kutta scheme (Appendix C). For case I, 

ALF{1) = 1/2, ALF(2) = 0, ALF(3) = 1 at the cold wall (x = 0) and 
ALF(4) - 1, ALF(5) = 0, ALF(6) = 2 at the hot wall (x = 1). Note 
that Y j to (equation (87)) is defined in computer language as Y(J), 
J = 1 to 6. If the values for Y (4), Y (5), and Y (6) at x = 1 are not 
equal to ALF(4), ALF(5), and ALF(6) within a specified error, 

0. 001, then an iterative scheme is devised to change the initially 
guessed values YO(4), YO(5), and YO(6) until | ALF(J) - Y(J)| < 

0. 001, J = 4, 5, 6. The details of the iterative scheme are given in 
Appendix C. 

The integration step size was changed to see if it had any 
effect on the result. Two schemes were tried: D = 0.001 across the 
entire flow field and D = 0. 001 near the two walls and D = 0. 01 in 
the rest of the flow field. The numerical results were essentially 
unchanged from those given by using D = 0. 01. 

Figures 3, 4, and 5 indicate how n 2 (0), u 2 (0), and |3 (0) change 
as Re increases from 0. 01 to 100 for case I. Similarly, Figures 6, 
7, and 8 show how 7^(1), u^l), and 8^1) behave over this Re range. 
It is worthwhile to note that n 2 (0), u 2 <0), and ^(0) approach limiting 
values 0. 6702, -.4105, and 1.553 respectively for Re >10. n (1) 
(Figure 6) is a bell shaped curve which has a peak value of O'. 83 5 at 
Re - 1.8 and approaches 0. 646 as Re - 100. IT (1 ) (Figure 7) has a 
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maximum value of 0. 405 at Re = 3.0 and decreases to 0. 37 at Re = 
100. Finally, (3^(1) (Figure 8) approaches 1.73 as Re -* 100. 

The behavior of the mean vapor velocity u(x) for case I is illus 
trated in Figure 9. At Re = 0.01 u is constant. As Re increases 
from Re = 0. 01 to Re = 1.0, the curve for u (x) has an essentially 
constant positive slope from x = 0 to x = 1. Physically this means 
the flow accelerated from the hot wall to the cold one. As the 
magnitude of Re becomes greater than one a point of inflection begins 
to appear in the u(x) curve. Finally, for Re > 6 the curve u(x) 
becomes concave upward except at the cold wall where it is concave 
downward. Physically, this corresponds to a flow which accelerates 
at the hot wall reaches a maximum value and then the decelerates 
toward the cold wall. At x = 0. 05 the flow starts to accelerate again, 
u is negative because the direction of integration (cold to hot) is 
opposite to that of the vapor motion. Velocity curves are plotted up 
to Re = 12.2. Figure 10 contains three of the u curves shown in 

i 

Figure 9 replotted on a scale comparable to that used in later 
figures. Figure 11 and Figure 12 indicate how p (x) and 8(x) behave 
for Re = 12.0. 

The curves for u^ (x) and u^ (x) are illustrated for Re = 11.7 
in Figure 13. Similarly n^x), n^fx) and g (x), 8 2 (x) are plotted in 
Figures 14 and 1 5 respectively for Re = 11.7. These figures 
indicate the existence of two regions of rapid change: one at the 
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cold wall and the other at the hot wall. 

Past Re = 12.0, it becomes increasingly difficult to use the 
method of integration discussed previously. However, at Re = 12.0, 
n^O), u^fO), and g^O) have almost reached their limiting values, i.e. , 
their values remain unchanged as Re -• 00 . This fact allows use of 
the integration technique of shooting-splitting. The forward marching 
scheme used is a fourth order Runge-Kutta with an integration step 
size D = 0. 0001 up to x = 0. 005 and D = 0. 001 past 0. 005. The 
integration process is begun at x = 0 for Re =' 100. 0 (Kn = 0. 0125 by 
(41)) and the values n^(0), u^(0), 8^(0), n^fO), 

Re = 12.0 are used to start the integration. As the integration 
moves forward u^ will eventually go to zero at some x < 1 and the 
integration process is stopped at that point. Next, the process is 
started over at x = 0 but for a value of u^fO) that is 0. 01 greater 
than u_(0) . . The new curve for u (x) eventually becomes 

w i\ 6 “ * 4 1 U 6 

greater than one at some x. u^ = 1 is used as a convenient cutoff 
point. If one takes the average between these two values for u^(0) 
and repeats the process either an up curve or a down curve will 
result. The averaging process is repeated until u (0) for the up 
curve and down curve match to five decimal points. Then at a point 
x = .15 (arbitrarily selected. Figure 16), the values of Y. for the 
two curves are averaged and a new curve is found. This curve 
either goes up or down, but more importantly, it extends further 


u^(0) and 8_(0) for 
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out into the flow field than any of the curves starting at x = 0. If this 

new curve is up for the sake of argument, then its values of Y and 

i 

those of the down curve which starts at x = 0 are averaged at x = 0. 15 
and a new up or down curve is calculated. In this manner the flow 
behavior for large Re can be found (Figures 16, 17, 18, 19, 20, 21, 
22, 23 , and 24). 

Each point at which ^(x) = 0 can be thought of as a position of 
the hot wall and the corresponding value of Re is found by multiplying 
x • 100 since Re is linear in distance. Furthermore, the shooting- 
splitting technique can be carried out pastx >1 and the resulting Re 
would be greater than 100. F rom Figures 18, 16, and 22 the various 
values of r^, u^ and ^ corresponding to 0 are plotted on 
Figuies 6, 7, and 8. Note that x = 1 corresponds to Re = 100 in the 
scale used in Figure 16 to Figure 24. However, each Re corres- 
ponding to u 2 (x) = o can also have x = 1 just by redefining the scale 
since it is linear. 

Figures 16 to 24 indicate that an equilibrium situation is 
attained at x = 0.2. Past this point u = u j= u^ n^ n, and 

8 - 2T as the limits in the above figures illustrate. As Re 
becomes large the regions of rapid change become thinner (Figures 9 
and 18). In the large Re limit the flow field has the following 
behavior: the velocity u accelerates from its values at the hot plate 

(x = 1) to an equilibrium speed and then decelerates from equilibrium 
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as the vapor approaches the cold wall. Very near this plate u 
accelerates again; both the density and temperature at x = 1 decrease 
to their respective values in the equilibrium state whereas near the 
condensing surface the temperature decreases and the density 
increases to a peak value and then decreases slightly. For Re -* 00 , 
the total enthalpy (Figure 2 5) decreases slightly at the hot wall and 
then approaches equilibrium. Near the condensing surface, it 
decreases strongly.' This flow pattern for large Re is like the result 
found by Collins and Edwards [ 16 ] except they found H to be constant 
at the evaporating surface. In the low Re limit (Re = 0.01) u", p", "p, 
etc. , are found to be essentially constant. 

In Figure 26 the velocity u is plotted for case III (92c). For 
Re = 0. 01 u is constant. In the range of Re between 0. 1 and 1,0 the 
vapor velocity u decelerates over most of the flow field and then 
accelerates very near the cold wall. This behavior is completely 
different from that of case I in the same Re range. For Re > 1.0 
the u curves follow the same pattern as those in Figure 12 (case I). 
Furthermore, n^O), u^(0), and jB^O) approach limiting values around 
Re - 14 and the shooting-splitting technique can be used. 

Case II (92b) has boundary conditions exactly opposite to those 
of case III. The mean velocity IT for case II (Figure 27) exhibits the 
same behavior as u in case I (Figure 9) except that there is no 
acceleration at the hot wall. 




Figure 
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Finally, from equation (86), there results 


Re - /"y — — • 

pu V7 rt7 


3 /RT d 


= / y M R( 


Y M Ke 
u 


where u - mean velocity. Since M ~ 0(1) and y = 5/3, Re ~ 0(Re). 

u pu 



Chapter 7 


CONCLUSIONS 

The two plate problem is solved for a monatomic vapor 
composed of Maxwell molecules. Lees' moment method is used to 
obtain a set of six non-linear moment equations whose solution, 
subject to the boundary conditions of this problem, is continuous 
over the range of flow conditions from free molecular to continuum. 

To obtain an analytical solution to this problem, small devia- 
tions from equilibrium are considered. A first order perturbation 
analysis is used and the mean velocity u is the small parameter in 
the problem. For Kn -* 00 , the flow properties and the vapor 
velocity are constant between the two plates and the value of the 
evaporation coefficient is one. When Kn ~*0, the flow is either 
density or temperature dominated depending on whether AT / AN < 
2.19 or AT/AN > 2. 19 respectively. The evaporation coefficient in 

this limit is 0.83. The Re based on the mass flux at the cold wall 

pu 

can be small or large depending on the magnitude of d/v. 

The two plate problem for strong non-equilibrium conditions 
is solved as an initial value problem. The direction of integration 
is opposite to that of the vapor motion. Three cases are considered. 
For Re -0 the vapor properties and velocity are constant across the 
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flow field. As Re -* 00 , two regions of rapid change appear: one at 
the hot wall and one at the cold wall. The flow accelerates at the hot 
wall to an equilibrium velocity, decelerates from equilibrium to a 
local minimum velocity as it approaches the cold wall, and then 
accelerates again. The vapor temperature and density decrease at 
the hot wall to their respective equilibrium values. However, at the 
cold wall the vapor density increases rapidly, reaches a peak and 
then drops slightly whereas the vapor temperature decreases. 
Finally, the total enthalpy decreases slightly at the hot surface to an 
equilibrium value and then drops sharply at the cold wall. 
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APPENDIX A 


DERIVATION OF EQUATIONS OF TRANSFER FOR Q(c.) = me. 

^ J 


As previously discussed, Maxwell's inverse fifth power law of 
molecular repulsion is used in evaluating the collision integral. For 
such molecules the collision integral becomes 


where 


AQ 


- /(mj+ m^K J j J dc dc 2 » 


“ ' 2tt 


J = J J (Q'-Q) de a da 


o o 


(A. 1 ) 


(A. 2) 


and a and e are geometric parameters describing the collision 
process. The description and interpretation of the parameters 
appearing in equation (A. 2) and the subsequent development have been 
treated by a number of authors. For details the reader is referred 
to the work of Lees [ 9 ] whose nomenclature and methodology is 
adopted in this work. 

For the moment of interest the difference in Q resulting from 
a collision is 


Q'-Q 

m 



(A. 3) 


To express Q'-Q in terms of c., use is made of the expression for c' 

J j 
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originally developed by Maxwell [9] , 


c . = 


2 / 2 2 

c.+ (c -c). cos (072) + V V - (c -c). sin 0' cos (e+ U) . ) 
J J 1 J jk 


(A. 4) 


where 0', e, and u> are angles in various planes which describe the. 
binary collision. By substituting this expression into (A. 3), Q'-Q 
becomes 


Q'-Q 2 2 2 2 2 

— — — - 3c. a + 3c. a 1 cos e 1 + 3c. a + 6c.aa'cos e 1 + 3c a' cos e 1 
m J J J J J 


3.2 


2 2 


+ a + 3a a' cose' + 3aa' cos s' + 


,3 3 , 

a 1 cos e 1 


(A. 5) 


where 


a = (c ^ - c). cos (072) 
J 


a = 


1 / 2 ,2 

2 V V " sin 01 

J 


e' = 


e + 

Jk 


After placing this expression for Q'-Q into equation (A. 2), integration 

over e is performed first. It is noted that terms proportional to odd 

powers of cos me or sin me (m ^ 0) integrate to zero over the 2 tt 

2 

limit hence terms involving cos e 1 and cos e' integrate to zero. 

2 

Using the fact that integration of cos (e + ) is tt and rearranging, 

j k 

Q'-Q 

the integral of over e becomes; 

m 
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(c'.^-c^)de = 2 tt -( (c -c).[" c 2 + c .c + c 2 1 cos 2 (0'/2) 

J J L 1 JL jl j jl j J 

- ? [ c j( I <c r c, j 2 • 1 y2 ) + (c i + c V c r c) j 1 sin e ' 

- ( C 1 ~ c )j [ \ ^ C 1~ C ^ “ F v2 J sin2 ®' cos 2 (6'/2) | . (A. 6a) 


Since this expression must be symmetrical with respect to the probe 
and colliding gas molecules, an equivalent statement of the bracketed 
term in equation (A. 6a) is 


{ \ = (c-c ).[ c 2 + c.c + c 2 ] cos 2 (6'/2) 

^ J 1 J Jl J Jl J 

•K C jl(l (c r C, j ' l y2 j + lc+cpjfcj-cl^.ln 2 *' 


- (c-Cj)^ — (Cj^-c) 2 - |- v 2 J sin 2 0' cos 2 (0'/2) . (A. 6b) 


To obtain symmetry in the expression for J the bracketed terms in 
equation (A. 6a) and (A. 6b) are summed and divided by two to yield 


T m , , 

J = 4 < c i + c 




3 V 2 - 2 (C 
2 ' 1 


CO 

^J" 


2 

sin 0’ 


a da 


(A. 7) 


Maxwell has evaluated the integral in equation (A. 7) and the result is 



95 


p 2 

- tt sin 0' a da = 1.3682 


(A. 8) 


Thus, the equation for J is 


J = A 2 (c 1 +c).[|v 2 -|(c 1 - c)2] 


(A. 9) 


To allow physical interpretation at a later point in the develop- 

2 

ment it is convenient to express V in terms of thermal velocities 
since the mean velocity of probe and colliding particles are identical 


2 . 2 


v = £ (c -c). = £ (c + c -2c c). , 

i=l i=l 


(A. 10) 


In making this substitution the collision term for a single component 
gas may be written as 


AQ = |mA 2 /2mK J J ff^ + c). [| V 2 - | (Cj-C) 2 ] dc dCj. (A. 11) 


Evaluation of the integral using c.= C.+ u. yields 

J J J 


AQ = |mn 2 A 2 /2mK { | [ C lj C i + C j c2+ 2u j (C l + C 2 )] 


IK^-ycW-]} 


(A. 12) 
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Again, calling upon the symmetry between probe and colliding 
particles, there results 


AQ 


= imn 2 A /2^K T I C.C 2 + 3u. C 2 - |c 3 - 9u. C 2 1 
2 2 L 2 j j 2 j J J J 


(A. 13) 


To simplify this expression use is made of the equation derived by 
Maxwell for viscosity based on a local full range Maxwellian 
velocity distribution 


P = 


kT 


^ A^ v/ 2mK 


(A. 14) 


With this result the collision term becomes 


AQ 


= p £ T C.C 2 /2 + u. C 2 4 C 3 - 3u. C 2 ] . (A. 15) 

P l J j 2 j J j J ' ' 


In t,erms of the shear stress and heat flux it may be rewritten 


AQ = 2. [ q.+ 3u.T. . - \ p C 3 1 . 
P L J J JJ 2 p j J 


(A. 16) 


The corresponding moment equation is 


4r fc 3 dcl+4~ fmf f c.c 3 dc*l = - [q.+ 3 u.t..- y p C 3 ] 

St L J J dx. L J 1 j J M L J J JJ 2 j J 


(A. 17) 



APPENDIX B 


INTEGRALS USED TO EVALUATE MOMENT EQUATIONS 


f e" §2/P dS = i /ITT 


+ er£ (7F)] 


(B.l) 


f 


%e-^ d5 


3 -c /p 
2 G 


(B.2) 


= |[. e 3 ]*[l + e ri ( 7 |-)]-Iece- c2/ e (B.3, 


r c f 3 -? 2 /e 

S e d§ 


3 , 2 -c /3 

- - (3 + c ) e 


• (B. 4) 


C .4 -? 2 /g 


5v?r, . .A c NT Pr_3,3_i -C 2 /P 


S c d? = |[ np 5 f[i +erf( 7 |-)j -Ip+lpc^ 


(B.5) 
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APPENDIX C 


COMPUTER PROGRAM 

The six moment equations that are solved in the subroutine 
DRVTV are: 

continuity, i = 1 


rz — dn 1 _ rz — dn 2 

( u i E WV T,x i> — + <VW — 

dx dx 


du i _ dU 2 n i d ^1 n ? dB ? 

+ n 1 E 1 — + « E + 1/2 — X — - !/2 = - X, — = 0 

dx J tt Pi dx J ttB 2 dx 

(C.l) 


x -momentum equation i = 2 


_ r ~ ' ^ dn / B i d u 

( [ V 2 + Ui ] E l +u 1 x iVVv 7 t '( ffi ii/7V 2 Vi E i) — 

dx 


dx 


n E dj3 _ __ / — dn 

+ — -r + ( Ce 2 /2 +u 2 3 e 2- vA'Oj 

dx dx 


. . ^2 x d U 2 n 2 E 2 d ^2 

+ ( 2 u 0 n 0 E 0 - 2n 0 X 0 J — ) —3 = 0 

' dx ■ * dx 


C 2 U 2 n 2“2 ~“ 2"*2 


(C. 2) 
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x-momentum flux, i = 3 


dx 

r 2 n du i 

+ [ 3n i u i E i +1 - 5n i B i E i +3n i u i7 - x i] — 

dx 


r 3 3 - \\ -I d8 l r- / '_3 N bz - -2 -| dn ' 

+ [ 2 n l “iV 2 "l V V X lJ 7T + [(“2 +1 - 5u 2 B 2>2-V -( B 2 + U 2 )X 2 ]^ 

dx dx 


B 2_. T dB 2 


r _ 2 3 i? 1 du ? 3 r 

+ [3n 2 u 2 E 2+? n 2 B 2 E 2 - ~ X^ 

dx " dx 


_ ^ 


2Re 


n{2„u 2 -2B 2 + i (n lSl E 1+ „ 2 6 2 E 2 )} =0 


(C.3) 


y-momentum flux, i = 4 


r ri dn i - - du i r __ 3 _ /p 7 -, d B, 

K ( “lV V ~ V + “ 1 S 1 E 1 — + [ n u lEl 4 - . y v Xj — 

dx dx dx 


+ 


[ B 2 ( u 


3 2 dn _ 

2 E 2-V ~ X 2 ) — + Vz E Z — ♦ [n 2 u 2 E 2 - |n 2 ^ * 2 


do 


2 . r— 


dx 


dx 


3 - bz,, -I dp 2 

J dx 


2R e 


n[2n u 2 - 2 B 2 + ± <■>, f>,E, + n^E.,)] = 0 


(C.4) 
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x- energy flux, i = 5 


[( I + 4u f e i + u i >i + /^C u i + T u i B i) x i]-r 

dx 


[ n i7T-< 6B i +4u f> x i H "i (z "iV u i 3 > e i]-^ + 


dx 


r s 5 _ —2\ — — / 

[“lG 3 l +4u l) E l + 4u l u lV — X l]— + 

dx 

[ ( f ®2 + 4 “2 \ + V ) V y?C “2 + i Vz ) X 2 ] -3 

dx 


r _ I $ 2 _ _2 _ 3 du 

[ ~ n 2 V V < 6 B 2 + 4 u 2 >X 2 + 4 n 2 ( 2 u 2 9 2 + U 2 )E 2 ] ~ + 

dx 


[" 2 (fV- 2 > 2 - 4 WK] d -p 4 


Ren^2B,- u (B, + n, p , E, + n 0 f3 0 E, 


3 L 3 -i - 2 ’ *‘n“! ' “ 2 ,J 2~2 


>] ■ 0 


(C. 5) 


flux of C in x direction, i = 6 
x 


[(t B i + 3u iV u i 4 ) E i + y^ ( ^ 3+ f 4 

dx 


[n 1 u 1 (4u 1 2 + 6B 1 )E 1+ 4 nj^- (oft 0^ ] ^ + 


du. 

dx 


[ 3 n i( T- +u i 2 ) E i + 3 n i u i VV x i]-^ + 
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r r 3 — 2 — 2 ~ - 4 

[Cl P 2 + 3U 2 B 2 +U 2 


X'/t ‘“M 

dx 


+ 


[ n 2 U 2 (4u 2 2+6p 2 )E 2 - 4n 2 V # (U 2 2+B 2 )X 2 ] 


du. 


dx 


+ 


3 


[ 3 n ri + u 2 ^) 

L 2 V 2 2 J 


3o dP. 


+ u ) E„ - 3n U 


2 ’ “2 ;^2 _ J “2“2V V X z]~^ +2Re n [ +2B 3 


- U [ B 2 + 4 (n l 3 l E l + n 2 e 2 E 2>] - | n i B lC U l E l + y V X l) 


-f n 2 P 2( U 2 E 2 


“ X 2 


)] = o 


(C.6) 


where 


E i ~ 1. + erf Q u^ / /jT^ 

= 1 - u 2 //lj) 


E. 


X = e 


X 2 


-h 2/ h 


n 2 ^ n l E l + n 2 E 2^ 


u = 


_ n 2 U 2 E 2 +n l U l E l + ( H l X l V if " n 


2 X 2 V TT 


0 


n l E l + n 2 E 2 


The complete computer program is listed on the following 


pages . 
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C MAIN PROGRAM - -INTEGRATION FROM COLD TO HOT WALL 
IMPLICIT REAL *8 (A-H.O-Z) 

REAL *8 K,L,M 

DIMENSION ALF(6),YO(6),W(6), B(3),K(6,4), Z(6), 

Y{6), DYO(3), DY1 (3), Y OSTO(6), Y 1ST 0(6), C(3, 3), 

L (6), M(6), DYA (3), F(6, 6), G(6) 

INTEGER CTR1.CTR2 
CALL $TIME$ (45) 

3 READ (5, 10, END = 9999)(ALF(I), 1=1, 6), (YO(J), 

J = 4, 6), RE 
10 FORMAT (6E12.8) 

A - 1 . 0/DSQRT (3. 141592DC) 

15 CTR2 = 0 
20 X = 0. 0 
CTR 1 =0 

C X IS THE INDEPENDENT VARIABLE RANGING FROM 0 

TO 1 

C Y (I) 1=1 TO 6 REPRESENTS THE VARIABLES N,U, AND BETA 

Y (1 ) = ALF(l) 

Y(2) = ALF{2) 

Y (3 ) = ALF(3) 

Y (4) = YO{4) 


Y{5) = YO{5) 
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Y<6) = YO(6) 

C THE EVAPORATION COEFFICIENT = EVCOEF; THE TOTAL 
ENTHALPY = TOTH 

CALL PHY SPR (Y , ALF, RHO, U , P, TEMP, TA UXX, 
QX,TOTH,B, EVCOEF, &800) 

WRITE (6,30) RE, (Y(I), 1=1 , 6), (B(I),I=1, 3) 

30 FORMAT (1H1 , 38X, 'FLAT PLATE EVAPORATION- 
CONDENSATION PROBLEM 1 
C/ / / 1H, 'INITIAL CONDITIONS RE = ', E12.4, ' 

Y (1) = ', E12.4, 

C'Y(2) = ', E12.4, ' Y (3 ) = ', E12.4, ' Y(4) = ', 

EI2.4 /1H, 

C21X, ' Y ( 5) = ', E12.4, ' Y(6) = ', EI2.4, ' 

B( 1 ) = ', E12.4, ' 

CB(2) = ', E12.4, ' B(3) = ', E12.4 lllm, 

'CTR1 X 

C D Y (1 )' , ' Y (2 ) 

Y(3) Y (4)' , 

C Y (5) Y (6)'/ 1H, 1 7X , 'RHO U 

CP ', 'TEMP TA UXX QX 

TOTH 

C EVCOEF' ) 

60 CALL DRVTV(Y, W,RE,F,G, &800) 
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C CTR 1 MULTIPLIED BY 100 GIVES THE POSITION X. 

C PRINT OUT OCCURS ONLY AT THE SPECIFIED POSITIONS 
GIVEN BY CTR 1 


IF 

(CTR1 

.EQ. 

0) 

GO TO 215 

IF 

(CTR 1 

.EQ. 

1) 

GO TO 215 

IF 

(CTR 1 

.EQ. 

2) 

GO TO 215 

IF 

(CTR 1 

.EQ. 

5) 

GO TO 215 

IF 

(CTR 1 

.EQ. 

10) 

GO TO 215 

IF 

(CTR 1 

.EQ. 

15) 

GO TO 215 

IF 

(CTR 1 

.EQ. 

20) 

GO TO 215 

IF 

(CTR1 

. EQ. 

2 5) 

GO TO 215 

IF 

(CTRi 

.EQ. 

30) 

GO TO 215 

IF 

(CTR1 

.EQ. 

40) 

GO TO 215 

IF 

(CTRI 

.EQ. 

50) 

GO TO 215 

IF 

(CTRI 

.EQ. 

60) 

GO TO 215 

IF 

(CTRI 

.EQ. 

70) 

GO TO 215 

IF 

(CTRI 

.EQ. 

75) 

GO TO 215 

IF 

(CTRI 

.EQ. 

80) 

GO TO 215 

IF 

(CTRI 

.EQ. 

8 5) 

GO TO 215 

IF 

(CTRI 

.EQ. 

90) 

GO TO 215 

IF 

(CTRI 

.EQ. 

95) 

GO TO 215 

IF 

(CTRI 

.EQ. 

98) 

GO TO 215 

IF 

(CTRI 

.EQ. 

99) 

GO TO 215 
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IF (CTR 1 .EQ. 100) GO TO 215 
IF (X .LT. 1.0) GO TO 217 

215 CALL PHYSPR (Y , A LF, RHO , U , P, TEMP, TA UXX, 

QX, TOTH, B, EVCOEF, &800) 

WRITE (6,216) CTRl.X.D, (Y (I), 1 = 1 , 6), RHO, U, P, 

TEMP, TAUXX , QX , TOTH, EV 
CCOEF 

216 FORMAT (1H, 14, 1P8E13.4/1H, 11X, 1P8E13.4) 

IF (X - 1.0) 217, 5070, 5060 

C INTEGRATE Y (J) TO X+D USING 4TK ORDER RUNGE-KUTTA. 

217 D = 0. 01 

C CORRECT LAST STEP FOR COMPUTER ROUND OFF ERRORS. 
IF ( (X+D) .LT. 1.0) GO TO 2 18 
D= 1 . 0 - X 
X=l. 0 
GO TO 220 

218 X=X +D 

CTR 1 = CTR 1 + 1 
220 CONTINUE 

NEXT EVALUATE THE COEFFICIENTS USED IN 
EXPRESSION FOR Y(J) 

DO 230 1=1,6 
K(I, 1) = D * W (I) 


C 
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230 Z(I) = Y(I) + 0. 50 * K(I, 1) 

CALL DRVTV(Z, W,RE,F,G, &800) 

DO 231 1=1,6 

K(I, 2 ) = D * W(I) 

231 Z(I) = Y (I) + 0. 50 * K(I, 2) 

CALL DRVTV(Z, W.RE.F.G, &800) 

DO 232 1=1,6 

K(I,3) = D * W{I) 

232 Z(I) = Y (I) + K(I,3) 

CALL DRVTV(Z,W,RE,F,G, &800) 

C NOW INTEGRATE Y(J) TO X+D USING THESE 

COEFFICIENTS 
DO 233 1=1,6 

233 Y (I) = Y (I) + (K(I, 1) + D * W(I) ) / 6.0D0 + 

(K(I, 2) + K(I, 3))/ 3.0D0 

IF (Y(3) .GE.0) GO TO 240 
WRITE (6,1) 

1 FORMAT (1HO, 10X, 'SQUARE ROOT OF NEGATIVE 
NUMBER Y (3 )’) 

WRITE (6.4) ((F (I, J), J=l, 6), W(I), G(I),I=1, 6) 

4 FORMAT (1X1P 6E 15.7) 

GO TO 800 


240 IF (Y (6) . GE. 0) GO TO 245 
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WRITE (6,2) 

2 FORMAT (1HO.10X, 'SQUARE ROOT OF NEGATIVE 
NUMBER Y (6 ) 1 ) 

WRITE (6,9) ( (F (I , J ) , J = l, 6), W(I), G(I),I=1 , 6) 

9 FORMAT (1X1P6E 15.7) 

GO TO 800 
24 5 CONTINUE 

IF (X-1.0) 60,60,5060 
5060 WRITE (6, 5062) X 

5062 FORMAT (1H1, 10X, 2 5 H X GREATER THA N 
1.0, X = , E 16. 7) 

5070 WRITE (6,5000) (B(I), 1=1,3) , (ALF(J), J=4,6) 

5080 FORMAT (1H, 'B(I), 1=1,3 = ' 1P3E13.4, 4X, 

'ALF(J), J=4, 6 = ', 1P3E - 
C13.4) 

C ITERATION TECHNIQUE: CHANGE YO(4) FIRST , 

C LEAVING YO(5) AND 70(6) UNCHANGED; USE THE 
C OLD VALUE OF YO(4), CHANGE YO(5), AND LEAVE 
C YO(6) UNTOUCHED; NEXT USE THE OLD VALUES OF 
C Y 0(4) AND YO(5) AND CHANGE YO(6). FINALLY 
C CHANGE YO(4), YO(5), AND YO(6) SIMULTANEOUSLY. 

C REPEAT PROCESS UNTIL ABS(ALF(I)-Y (I)) < ERROR. 



DO 6000 J=4,6 

IF {(DABS(ALF(J) - Y(J))) .GT. 0. 001) GO TO 
6010 

6000 CONTINUE 
GO TO 3 

6010 IF (CTR2 .GT. 0) GO TO 6020 
CTR2 = CTR2 + 1 

DO 6011 J=4, 6 
YOSTO(J) = YO(J) 

6011 YISTO(J) = Y (J) 

DYO(l) = (ALF (4) - Y (4 ) ) / 10.0 
DYO(2 ) = (ALF (5) - Y{5)) / 10.0 
DYO(3 ) = (ALF{6) - Y(6)) / 10.0 
DYA(l) = ALF(4) - Y(4) 

DY A(2) = ALF(5) - Y(5) 

DYA(3) = ALF(6) - Y{6) 

6012 YO(4) = YOSTO(4) + DYO(l) 

GO TO 20 

6020 IF (CTR2 .GT. 1) GO TO 6030 
CTR2 = CTR2 + 1 
DY1(1) = Y (4) - Y lSTO(4) 

DY 1 (2) = Y (5) - Y 1ST 0(5) 

DY1(3) = Y(6) - YlSTO(6) 

DO 6021 J=1 , 3 
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6021 C(J, 1) = (DY1(J) / DYO(l) ) 

YO{4) = YOSTO(4) 

6013 YO(5) = YOSTO(5) + DYO(2) 

GO TO 20 

6030 IF (CTR2 .GT. 2) GO TO 6040 
CTR2 = CTR2 + 1 

DY 1(1) = Y(4) - Y lSTO(4) • 

DY 1(2) = Y (5) - YlSTO{5) 

DY 1 (3 ) = Y(6) - Y lSTO{6) 

DO 6031 J=l, 3 

6031 C(J, 2) = (DY1(J) / DY 0(2 ) ) 

YO(5) = Y OSTO(5) 

6014 YO(6) = YOSTO(6) +DYO(3) 

GO TO 20 

6040 DY 1(1) = Y (4) - YlSTO(4) 

DY 1 (2) = Y(5) - Y lSTO(5) 

D Y1 (3 ) = Y (6) - Y lSTO(6) 

DO 6041 J=1 , 3 

6041 C(J,3) = (DY1{J) / DYO (3) ) 

C SUBROUTINE DMINV IS AN IBM MATRIC INVERSION 
C PROGRAM - DOUBLE PRECISION. 

CALL DMINV (C,3,Q,L f M) 

IF (Q . EQ. 0.0) GO TO 6070 



Y 0(4) = Y OST 0(4) + C(l, 1) * DYA(l) + C(l,2) * 
DYA(2) + C(l,3) * DYA (3) 

YO(5) = YOSTO(5) + C(2, 1) * DYA(l) + C(2,2) * 
DYA(2) + C (2 , 3 ) * DYA{3 ) 

YO(6) = YOSTO(6) + C(3, 1) * DYA(l) + C(3,2) * 
DYA(2) + C(3, 3) * DYA (3 ) 

GO TO 15 

6070 WRITE (6,6071) 

6071 FORMAT (1H1, 10X, 2 9HDETERMINANT OF C 
EQUALS ZERO) 

GO TO 9999 
9999 CALL EXIT 

C THE FOLLOWING INCREMENTAL CHANGES ARE MADE 
C ONLY IF Y(3) OR Y(6) BECOME NEGATIVE. THIS IS 
C DONE IN ORDER TO RESTART THE PROGRAM AGAIN. 

800 IF (CTR2 . EQ. 0) GOTO 9999 

IF (CTR2-2) 801, 802, 803 

801 DYO(l) = -DYO(l)/2. 

ITEK = ITEK + 1 

IF (ITEK .EQ. 5) GO TO 9999 
GO TO 6012 

802 DYO(2) = -DYO(2)/2. 


ITEK = ITEK + 1 



IF (ITEK. EQ. 5) GO TO 9999 
GO TO 6013 

803 DYO(3) = -DYO(3)/2. 

ITEK = ITEK + 1 

IF (ITEK. EQ. 5) GO TO 9999 

GO TO 6014 

END 

C SUBROUTINE PHYSPR 

SUBROUTINE PHYSPR (Y , ALF, RHO, U, P.TEMP, 
TAUXX, QX, TOTH, B, EVCOEF, *) 

C PURPOSE - TO EVALUATE PHYSICAL PROPER- 

TIES IN THE FLOW FIELD 
IMPLICIT REAL *8 (A-H.O-Z) 

DIMENSION Y(6),ALF{6), B(3) 

A = 1.0 /DSQRT (3.141592D0) 

IF (Y (3). GE. 0) GOTO 1000 
WRITE (6, 1001) 

1001 FORMAT (1H0, 10X, 'SQUARE ROOT OF A 
NEGATIVE NUMBER Y(3) SUB P') 

RETURN 1 

1000 IF(Y(6).GE. 0) GO TO 1002 
WRITE (6, 1003) 

1003 FORMAT (1H0, 10X, 'SQUARE ROOT OF A 



NEGATIVE NUMBER Y(6) SUB P') 
RETURN 1 

1002 R 1 = DSQRT (Y(3)) 

R2 = DSQRT <Y(6)) 

El = 1.0 + DERF (Y (2) / Rl) 

E2 = 1.0 + DERF (- Y (5) /R2) 

XI = DEXP (- (Y (2 ) ** 2)/Y(3)) 

X2 = DEXP (- (Y (5) ** 2) /Y(6)) 

ARX 1 = A * Rl * XI 
ARX2 = A * R2 * X2 
Y2S = Y(2)**2 
Y5S = Y(5)**2 

RHO = 0. 50 * (Y (1 ) * El + Y (4) * E2) 

IF (RHO. LT. 0.0) RETURN 1 

B(l) = . 50 * (Y(1)*(Y(2)*E1+ARX1) +Y(4) 

(Y (5) * E2-ARX2)) 

U = B(l) / RHO 
US = U**2 

UMY2S = (U - Y(2))**2 
UMY 5S = (U - Y(5))**2 
UX2MY2 = 2. 0 * U - Y(2) 

UX2MY5 = 2. 0 * U - Y(5) 

P = (Y (l)/6. 0D0) * ((UMY2S + 1. 50^Y (3 )) 



El - ARX1 * UX2MY2 ) 

C+ (Y(4)/6. ODO) * ((UMY5S + 1 . 50*Y (6))*E2 + 
ARX2 * UX2MY5 ) 

TEMP = P/RHO 

TAUXX = -(1.0D0/3.0D0)*(Y(1)*(UMY2S*E1 - 
UX2MY2 * ARX1 ) + Y(4) * 

C(UMY5S*E2 + UX2MY5 * ARX2 ) ) 

QX = .2 50 * (Y(l)*(-2. 50*Y(3) + UMY2S ) * 
(U-Y(2))*E1 + (2. 0*Y (3) + 

CY2S + 3.0* (US - U*Y(2)) )*ARX1) - Y(4)* 

((2. 50 * Y (6) + UMY5S ) * 

C(U-Y(5))*E2 4 (2. 0*Y (6) + Y 5S +3.0 -(US - 
U*Y(5)) )* ARX2 ) ) 

TOTH = 2. 50 * TEMP + . 50 * US 

SIGMXX = TAUXX - P 

B(2) = RHO-US - SIGMXX 

B(3 ) = RHO*U * (1. 50 * TEMP + . 50* US) - 

U * SIGMXX + QX 

EVCOEF = 2. 0 * B(l) / (A* (ALF(l) * (DSQRT 
(ALF(3))) - ALF (4) * (D 
CSQRT(ALF(6))) ) ) 

RETURN 


END 



C SUBROUTINE DRVTV 

SUBROUTINE DRVTV (Y , W , RE, F, G , *) 
IMPLICIT REAL *8 (A-H.O-Z) 

C PURPOSE 

C OBTAIN VALUES OF DERIVATIVES WJ AT 

GIVEN X KNOWING Y J AND RE 
DIMENSION Y(6), W(6), B (3 ), G (6), F (6, 6) 

A = 1. 0/DSQRT(3. 141592D0) 

IF (Y(3) . GE. 0) GOTO 1000 
WRITE (6, 1001) 

1001 FORMAT (1H0.10X, 'SQUARE ROOT OF A 
NEGATIVE NUMBER Y(3) SUB D 1 ) 

RETURN 1 

1000 IF (Y (6) . GE. 0) GO TO 1002 
WRITE (6, 1003) 

1003 FORMAT (1H0, 10X, 'SQUAREROOT OF A 
NEGATIVE NUMBER Y (6) SUB D') 

RETURN 1 

1002 R 1 = DSQRT (Y(3)) 

R2 = DSQRT (Y (6)) 

El = 1.0 + DERF (Y (2 )/R 1 ) 

E2 = 1.0 + DERF (- Y (5)/R2 ) 

XI = DEXP (- (Y (2) **2) / Y (3 ) .) 


X2 = DEXP (- (Y (5) **Z) / Y(6) ) 
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ARX1 = A * R 1 * XI 

ARX2 = A * R2 * X2 

YAXDR1 = Y(l) * A * XI / R1 

YAXDR2 = Y (4) * A * X2 / R2 

Y IE = Y (1 ) * El 

Y4E = Y (4) * E2 

Y2E = Y (2) * El 

Y5E = Y (5) * E2 

Y 12 = Y (1) * Y (2) 

Y45 = Y (4) * Y (5) 

Y2S = Y(2)**2 
Y5S = Y(5)**2 
RHO = . 50 * (Y1E + Y4E) 

B ( 1 ) = 0. 5 * (Y (1) * (Y2E + ARX 1 ) + Y(4) * 

(Y5E - ARX2) ) 

U = B (1 )/RHO 
USQ = U**2 

B (2 ) = . 50 * (Y (1) * ( { .50 * Y (3 ) +Y2S) * 

El + Y (2) * ARX 1 ) + 

C Y (4) * ( ( . 50 * Y(6) + Y5S ) * E2 - Y(5) * ARX2)) 
B{3) = .2 50 * (Y (1) * ( ( 2. 50 * Y(3) + Y2S ) * 

Y2E + ( 2. 0 * Y(3)+ Y2S ) * 

C ARX 1 ) + Y(4) * ( ( 2. 50 * Y(6) + Y5S) * Y 5E - 


C ( 2.0 * Y (6) + Y5S) * ARX2 ) ) 

F(1 , 1) = Y2E + ARX1 

F ( 1 , 2 ) = Y1E 

F(l, 3) = 0. 5 * YAXDR1 

F(1 ,4) = Y5E - ARX2 

F(1 , 5) = Y4E 

F(1 , 6) = -0. 5 * YAXDR2 

F (2 , 1) = (. 50 * Y (3) + Y2S) * El + Y(2) * ARX1 
F (2 , 2 ) = 2.0 * Y (1) * {Y2E + ARX1) 

F (2 , 3 ) = . 50 * Y1E 

F(2,4) = (. 50 * Y (6) + Y5S) * E2 - Y(5) * ARX2 
F(2, 5) = 2.0 * Y (4) * (Y 5E - ARX2 ) 

F (2 , 6) = . 50 * Y4E 

F (3 , 1) = (Y2S + 1.5* Y(3))*Y2E + ARX1* 

(Y (3 ) + Y2S) 

F (3 , 2 ) = 3.0 * Y (1) * ((Y2S + 0. 50 * Y(3)) * 

El + Y (2 ) * A RX 1 ) 

F(3,3) = 1. 5 * Y(l) * (Y2E + ARX1) 

F(3 , 4) = (Y 5S +1.5* Y (6)) * Y5E - ARX2 * 

(Y (6) + Y 5S) 

F(3, 5) = 3.0 * Y (4) * ({Y5S + 0. 50 * Y(6))* 

E2 - Y (5) *ARX2) 

F(3 , 6) = 1. 5 * Y (4) * (Y5E - ARX2) 
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F (4 , 1) = Y (3 ) * (Y2E + ARX1) 

F (4 , 2 ) = Y (3 ) * Y IE 
F (4 , 3 ) = Y(l) * (Y2E + 1. 5 * ARX1) 

F(4, 4) = Y (6) * (Y5E - ARX2 ) 

F (4 , 5) = Y (6) * Y4E 

F(4,6) = Y (4) * (Y5E - 1. 5 * ARX2 ) 

F(5, 1) = ( (1.25 * Y(3) + 4.0 * Y2S ) * 

Y(3) + Y 2S**2 ) * El 
C + Y(2) * ARX 1 * (Y2S + 3. 50 * Y(3) ) 

F (5, 2 ) = Y (1 ) * ARX 1 * ( 6. 0 * Y(3) + 4.0 * Y2S) 
+ 4. 0 * Y1E * Y (2 ) * (2. 0 * Y(3) + Y2S ) 

F(5, 3 ) = Y IE * (2.5* Y (3 ) +4.0 * Y2S) + 

4. 0 * Y12 * ARX 1 

F(5, 4) = (Y (6) * (1.25 * Y(6) + 4. 0 * Y5S) + 
Y5S;**2 ) * E2 - ARX2 
C * Y (5) * (Y5S +3. 50 * Y(6) ) 

F(5, 5) = -Y (4) * ARX2 * (6. 0 * Y(6) + 4. 0 * 

Y5S) + 4. 0 * Y4E 
C * Y (5) * (2. 0 * Y(6) + Y5S ) 

F(5, 6) = Y4E * (2. 5 * Y(6) + 4. 0 * Y5S) - 4. 0 * 
Y4 5 * ARX2 

F (6 , 1) = (Y (3 ) * (0. 75 * Y(3) +3.0 * Y2S) + 
Y2S**2 ) * El + 



C Y(2) * ARX1 * ( Y2S +2.5* Y(3) ) 

F(6, 2) = Y 12 * (4. 0*Y2S + 6. 0*Y(3) ) * El + 

4. 0 * Y (1 ) * ARX1 * {Y2S + Y (3 ) ) 

F(6, 3) = 3.0 * Y IE * {0. 5*Y(3) + Y2S) + 

3.0 * Y 1 2 * ARX1 

F(6 , 4) = ( (.7 50 * Y (6) + 3.0 * Y 5S) * Y(6) + 

Y 5S**2 ) * E2 - Y (5) 

C * ARX2 * ( Y5S + 2. 5 * Y(6) ) 

F (6 , 5) = Y4 5 * (4. 0*Y5S + 6. 0 * Y(6) ) * E2 - 
4. 0 * Y (4) * ARX2 * (Y 5S + Y(6) ) 

F(6, 6) = 3.0 * Y4E * (0. 5 * Y(6) + Y5S) - 
3.0 * Y45 * ARX2 
G (1 ) =0.0 
G (2) = 0. 0 

G (3 ) = -(2. 0D0/3. 0D0) * RHO * RE * (2.0 * 
RHO * USQ - 2.0* 

CB (2 ) + . 50 * (Y IE * Y{3) + Y4E + Y(6) ) ) 

G (4) = -G (3 ) 

G(5) = (4. 0D0/3. 0D0) * RHO * RE * (2.0 * B(3) 
- U * (B(2) + Y IE 
C* Y (3 ) + Y4E * Y (6) ) ) 

G (6) =2.0* RHO * RE * (-U *.(B(2) + . 2 50 * 


(Y1E * Y (3 ) + Y4E 
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C * Y(6) ) ) + 2.0 * B (3 ) - . 750 * Y(l) * Y(3) * 
(Y2E + ARX 1 ) - 

C. 7 50 * Y (4) * Y (6) * <Y5E - ARX2 ) ) 

C SUBROUTINE DSIMQ IS AN IBM PROGRAM WHICH 
C .SOLVES A SYSTEM OF SIMULTANEOUS LINEAR 
C EQUATIONS - DOUBLE PRECISION. 

49 CALL DSIMQ (F,G, 6, KS) 

DO 50 1=1, 6 

50 W(I) = -G (I) 

IF (KS. EQ. 0) GO TO 51 
WRITE (6, 52) 

52 FORMAT (1HO, 10X, 'KS IS ONE SINGULAR 
SOLUTION') 

RETURN 1 

51 CONTINUE 
RETURN 
END 

C SHOOTING-SPLITTING MAIN PROGRAM WHICH 
STARTS AT X=0. 

IMPLICIT REAL *8 (A-H.O-Z) 

REAL *8 K, L,M 

DIMENSION ALF(6), YO(6), W(6), B(3),F(6, 6), 

G(6),Y(6),Z(6),K(6,4) 
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INTEGER CTR1 


READ (5,10) (ALF(I),I = 1, 6), YO(4), YO(6),RE 


10 

3 

31 


FORMAT (6E12.8/3E12.8) 
READ (5,31, END = 9999) YO(5) 
FORMAT (E12.8) 


^i(o) is the critical 
in shooting-splitting 
It is constant near x= 
drops off sharply at > 


A = 1. 0/DSQRT(3. 141 592D0) 


20 X=0.0 


CTR1 =0 

Y (1 ) = ALF(l) 

Y (2) = ALF (2) 

Y (3 ) = ALF (3) 

Y (4) = Y O (4 ) 

Y (5) = YO(5) 

Y (6) = Y 0(6) 

CALL PHYSPR (Y , A LF , RHO, U, P, TEMP, 

TAUXX, QX, TOTH, B , EVCOEF, &800) 

WRITE (6,3) RE, (Y (I),I=1 , 6), (B(I), 1=1 , 3) 

30 FORMAT (1H1,38X, 'FLAT PLATE EVAPORA TION- 
C ONDENSATION PROBLEM 1 
C///1H, 'INITIAL CONDITIONS RE = ', 

E 1 2 . 6, ' Y (1 ) = ', E12.6, 


C' Y (2) = ', E12.6, ' Y (3 ) = ', E12.6, ' 


parameter 
technique 
0 and 
= 1 . 


Y (4 ) = ', E12. 6 /H , 
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C21X, * Y (5) = \ E12.6, 1 Y(6) = E12. 6 , ' 

B(1 ) = E12. 6, 

CB (2) = \ E12.6 , ' B(3) = E12.6 ///1H, 

'CTRl X 

C D Y (1 )' , 1 Y (2) Y (3) Y (4)' , 1 

C Y (5) Y(6)'/1H, 1 7X, 'RHO U 

CP ', 'TEMP TAUXX QX TOTH', 1 
C EVCOEF 1 ) 

60 CALL DRVTV(Y, W,RE,F,G, &800) 

IF (CTRl .EQ. 10000) GOTO800 

215 CALL PHYSPR(Y, ALF, RHO, U, P, TEMP, 

TAUXX , QX, T OTH, B, EVCOEF, &800) 

WRITE (6,216) CTR 1 , X, D, (Y (I), 1=1 , 6), 

RHO, U, P, TEMP, TAUXX, QX, TOTH, EVCOEF 

216 FORMAT (1H, 14, 1P8E13 . 4/ 1H, 1 IX, 1P8E13.4) 
INTEGRATE Y(J) TO X+D USING 4TH ORDER 
R-K AND VARIABLE STEP SIZE 

217 IF (CTRl .LT. 50) GO TO 81 
D = 0. 001 

GO TO 82 

81 D = 0. 0001 

82 X = X+D 


IF (CTRl .LT 


50) GO TO 91 



CTR 1 = CTR1+10 


GO TO 220 
91 CTR 1 = CTR1 + 1 

C NEXT EVALUATE THE COEFFICIENTS USED IN 

EXPRESSION FOR Y(J) 

DO 230 1=1,6 
K(I, 1) = D * W(I) 

230 Z(I) = Y(I) + 0. 50 * K(I, 1) 

CALL DR VT V (Z, W,RE,F,G, &800) 

DO 231 1=1,6 
K(I,2) = D * W(I) 

231 Z(I) = Y(I) + 0. 50 * K(I,2) 

CALL DRVTV(Z, W.RE.F.G, &800) 

DO 232 1=1,6 
K(I, 3 ) = D * W(I) 

232 Z(I) = Y(I) + K{I,3) 

CALL DRVT V (Z, W,RE,F,G, &800) 

C NOW INTEGRATE Y{J) TO X+D USING THESE 

COEFFICIENTS 
DO 233 1=1, 6 

233 Y (I) = Y (I) + (K(I, 1 ) +D*W(I) )/6. 0D0 + (K(I,2) 

+ K(I, 3 )) / 3 . 0D0 
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IF (Y(5) .GE. 0.0) GO TO 61 
IF (Y (5) . I.E. -1.0) GO TO 62 
IF (Y(3) .GE. 0) GO TO 240 
WRITE (6, 1) 

1 FORMAT (1H0, 10X, ’SQUARE ROOT OF 
NEGATIVE NUMBER Y (3 ) ' ) 

GO TO 800 

240 IF (Y(6) .GE. 0) GO TO 245 
WRITE (6,2) 

2 FORMAT (1H0, 10X, 'SQUARE ROOT OF 
NEGATIVE NUMBER Y (6 ) ') 

GO TO 800 
245 CONTINUE 

IF (X-1.0) 60,60,800 
9999 CALL EXIT 
800 WRITE (6,69) (B(I),I=1,3) 

69 FORMAT (1H0, 'B(I), 1=1,3 = ' 1P3E15. 6) 

GO TO 3 

61 WRITE (6,63) 

63 FORMAT (1H0, 10X, 'Y (5) IS GREATER THAN 
OR EQUAL TO ZERO ') 

WRITE (6,65) (B(I), 1=1 , 3) 

65 FORMAT (1H0, 'B(I), 1 = 1,3 = ' 1P3E16.6) 
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GO TO 3 

62 WRITE (6,64) 

64 FORMAT (1H0, 10X, 'Y(5) IS LESS THAN OR 
EQUAL TO MINUS ONE') 

WRITE (6,67) ( B (I ) , 1=1,3) 

67 FORMAT (1H0, ’B(I), 1=1 , 3 = ' 1P3E16.6) 

GO TO 3 
END 

C SHOOTING -SPLIT TING MAIN PROGRAM WHICH STARTS 
AT X = NUMBER. 

IMPLICIT REAL *8 (A-H.O-Z) 

REAL *8 K, L,M 

DIMENSION ALF(6), YO(6), W(6), B(3 ), F(6, 6), G (6),Y (6) ,Z(6),K(6,4) 
INTEGER CTR I 

C ALF(X) = VALUES OF UP CURVE AT X; Y0(J) = VALUES OF DOWN CURVE AT X. 

3 READ (5, 10, END = 9999) (A LF (I), 1=1 , 6), (YO(l), 1=1 , 6), 

CTR 1 , X, R E 

10 FORMAT (6E12. 8/6E12. 8/15, 2E12. 8) 

A = 1 . 0/DSQRT (3. 141592D0) 

DO 99 1=1,6 

99 Y (I) = (ALF(I) + YO(I))/2. 0D0 

CALL PHYSPR (Y , ALF, RHO, U, P, TEMP, 

TAUXX,QX,TOTH, B.EVCOEF, &800) 
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WRITE (6,30) RE, (Y (I) , 1=1 , 6), (B (I), 1 = 1 , 3 ) 

30 FORMAT (1H1, 38X, 'FLAT PLATE EVAPORATION 
-CONDENSATION PROBLEM 1 
C///1H, 'INITIAL CONDITIONS RE = ', 

E12. 6, ' Y (1 ) = ', E12. 6 , 

C' Y (2) = ', E12.6, * Y (3 ) = ', E12. 6, ' 

Y (4) = E12. 6 /1H , 

C21X, 'Y (5) = ', E12. 6, ' Y(6) = E12.6, ' 

B (1 ) = ', E12. 6 , 1 

CB(2) = E12.6 , ' B{3) = ', E12.6///1H, 

'CTR 1 X 

C D Y (1 )' , ' Y (2) Y (3) Y (4) 1 , ' 

C Y (5) Y(6)'/1H, 17X, 'RHO U 

CP 'TEMP TAUXX QX TOTH',' 

C EVCOEF' ) 

60 CALL DRVTV(Y, W, RE, F,G, &800) 

IF (CTR 1 .EQ. 20000) GO TO 800 

215 CALL PHYSPR(Y, ALF, RHO, U, P, TEMP, 

TAUXX, QX, TOTH, B, EVCOEF, &800) 

WRITE (6,216) CTR 1 , X, D, (Y(I), 1=1, 6), RHO, 

U, P, TEMP, TAUXX, QX, TOTH, EVCOEF 

216 FORMAT (1H, 15, 1 P8E1 3 . 4/ 1H, 1 IX, 1 P8E13 . 4) 
INTEGRATE Y(J) TO X+D USING 4TH ORDER 
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R-K AND VARIABLE STEP SIZE 
D = 0.001 

X = X+D 

CTR I = CTRI + 10 

C NEXT EVALUATE THE COEFFICIENTS USED IN 

EXPRESSION FOR Y(J) 

DO 230 1=1 ,6 

K ( 1 , 1 ) = D * W(I) 

230 Z (I) = Y (1) + 0.50 * K(1 , 1 ) 

CALL DRVTV(Z,W,RE,F,G,&800) 

DO 231 1=1 ,6 

K(l,2) = D * W(I) 

231 Z(I) = Y (I) + 0.50 * K(I,2) 

CALL DRVTV(Z,W,RE,F,G,&800) 

DO 232 1=1,6 


K(I,3) = D * W(I) 
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232 Z(I) = Y{I) + K(I,3) 

CALL DRVTV(Z, W.RE.F.G, &800) 

C NOW INTEGRATE Y(J) TO X+D USING THESE 

COEFFICIENTS 
DO 233 1=1,6 

233 Y (I) = Y(I) + (K(I, 1) + D*W(I))/6. 0D0 + 

(K{I, 2) + K(I, 3 ))/3 . 0D0 

IF (Y (5) .GE. 0. 0} GO TO 61 
IF (Y (5) . LE. -1.0) GO TO 62 
IF (Y (3) .GE. 0) GO TO 240 
WRITE (6, 1) 

1 FORMAT (1H0, 10X, 'SQUARE ROOT OF 
NEGATIVE NUMBER Y{3)') 

GO TO 800 

240 IF (Y(6) .GE. 0) GO TO 245 
WRITE (6,2) 

2 FORMAT (1H0, 10X, 'SQUARE ROOT OF 
NEGATIVE NUMBER Y(6) ') 

GO TO 800 
245 CONTINUE 

IF (X-2.0) 60,60,800 
9999 CALL EXIT 
800 WRITE (6,69) (B(I), 1 = 1,3) 



69 FORMAT (1H0, 'B(I), 1=1,3 = ' 1P3E15. 6) 

GO TO 3 

61 WRITE (6,63) 

63 FORMAT (1H0, 10X, 'Y (5) IS GREATER THAN 
OR EQUAL TO ZERO') 

WRITE (6,65) (B(I), 1=1,3) 

65 FORMAT (1H0, 'B(I), 1=1,3 = ' 1P3E16.6) 

GO TO 3 

62 WRITE (6,64) 

64 FORMAT (1H0, 10X, 'Y(5) IS LESS THAN OR 
EQUAL TO MINUS ONE') 

WRITE (6,67) (B(I), 1 = 1 , 3 ) 

67 FORMAT (1H0, 'B(I), 1=1,3 = 1 1P3E16.6) 

GO TO 3 


END 



APPENDIX D 


A NAVIER-STOKES TYPE FORMULATION OF THE 
LEES' MOMENT EQUATIONS 


The six moment equations to be used are 


continuity 


t( u i v/5 


(33) 


x-m omentum 


7 (n I B l E l +n 2 B 2 E 2> + KVJWI^ 


2 / B 2 n 

+ n., u n E„ - n„ u ./ — ) = B„ 

tt y 2 : 


2 22 " 2 2^2 


(34) 


energy 


|( n 1 a 1 3 Ei+n 2 . 2 


u^E,) + | (n x UjBjE^ n 


2 u 2 P 2 E 2 ^ 


+ 



n 2 U 2 2x 2 



) + K ; 


2 e i n i x i 




(3 5) 
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x-momentum flux 


!r[ n i (u i 3+ l B i u i ,E i +n i 


B 


>i +u i 2 > x i 


+ n 2 (U 2 + lh a 2 )E 2 n 2j-^ < B 2 +u 2 )X 2 ] 


+ I u [' 2n “ 2 + 2 V I (n l S l E l + n 2 B 2 E 2>] = ° 


(36) 


x-energy flux 


d 

dx 


["i (I b i 2 + 4 ii i 2 b i + “i B !> x i 


+ n z(l B 2 2+ 4u 2*2 + u 2 T (u 2 3+ I U 2 p 2 ) X 2 ] 


+ f ?[ 2B 3- U V “'VlVWV] = ° 


(38) 


x-flux of C 


x 


d^[ n l (? B l + 3B l"l + “1 >1 + VT (U 1 + 2 B 1“1> X 1 

+ n 2C 4 B 2 + 3B 2 U 2 + U 2 ) E 2 _n 2 JlT ^2 + 2 B 2 U 2* X 2] 
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+ 2 m[ - U(B 2 + ? (n i e i E l +n 2 B 2 E 2>> +2B : 


- 4 ‘“iVlV^W 


(Vi x i fi -Wz /¥) 


= 0 


and 


n = 


2< n lV n 2 E 2> 


T = 


-L [ t{ ( (u - u i >2+ iW /r < 2 u - u i) x J 


+ -r { ( (U -V 2+ 1 WV t- < 2u - V x 2 } ] 


- U | 2/ B, - U 2 /P ? 

where R is the gas constant, X = e ,X=e , E = 

1 2 1 

1 + erf / /p ~ , and E^ = 1 + e rf (- / /"j^j ) . These equations 

be rewritten in terms of speed ratio S = — . The result is 


/2RT 


continuity 


n i U l 


x 1 X 1 ^ , n 2 U 2^_ 1 “2\ 

^ 1 /s s i ^ 2 ^ 2 " ^ B ‘ 


X, 


x-m omen turn 

2 


n i u i 


Ej 2X 

<-2 +2E l + — , 

Sj S /tt 


, n u E 2X 

^♦-TrC^-v^O 

S 2 V" 


= B, 


(3 9) 


(D.l) 


(D.2) 


can 


(D.3) 


(D.4) 
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energy 


n l U l 


- E X, 2X 

C E l + 2 ~J + + — ) 

v 1 2 s 2 s ^ s 3 /^ y 


n 2 U 2 


CvlJ 

b 2 


X„ 2X„ 

2 2 N 

s Jvi s 3 /^ y 


x-m omen turn flux 




, _ E X. X. 
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b) 
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1 f 111, 222 
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)] = 0 


x-energy flux 


S 1 S 1 S 1 S l /n 


, ,, 4 z' 5 1 , 4 , 4^. 7 1 N 2 1 

n 2 2C4 4 + 2 + 0 E 2 " n 2 U 2 C 2 2) /- ] 

S 2 S 2 S ~ S - ^ 


2 2 


+ i £[2 

3 u L 


2 2 

„ r U i n i E i , u 2 n 2 E 2^-] 

uB 2" u ( —2 + — )J 

S 1 S 2 
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x-flux of C 


df 4/3 1 3 4/ 5 1 n X i 

dx L 1 1 V 4 4 2 J l 1 1 V 2 2A r 


S ! S 1 
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X, 
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(D. 8) 


and 


T = 


J_ f 

- n i u f E i 

2 

n l U ! 

( ? JL_ A 

nR l 


6 S i/n 

l A ' 


, n 2 U 2 2E 2 

2 
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( 2 u i "1 

A 

6s 2 a; 

^ u 2 • V X 2j 


(D.9) 


To develop a Navier-Stokes type representation from equations 
(D.3) to (D. 9), a procedure analogous to that used by Chapman and 
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Enskog is followed here. Let 

n x = n(l + n'p 
n 2 = n(l + n*) 
s x = S(1 + S*) 
s 2 =s(i+sp 
= u(l + Up 

u 2 U(1 + U 2^ 

where the differences 

n r n 2 = n(n i ■ n 2> 
u r "2 = u(u i ■ u 2> 
S 1 - S 2 = S < S I - S I> 


express a small deviation from equilibrium. By placing these 
expansions into equations (D. 1) and (D.3) to (D. 9), there results 

continuity 

nu = B^ (D. 10) 


x -momentum 


nu 
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S/W L 1 





(D. 11) 
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energy 


-S 2 c 2 

nu / 5 s 1 1 n. o v <•' "* ^ 
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(D. 12) 


x-m omentum flux 
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x- energy flux 
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(D. 14) 


x-flux of C 


x 


d . 6 3 


d~C 2 + rz + .4 


-S 


^)nu 4 = -2 — nu 3 • - — 
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(D. 15) 
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where the combinations 


density 


n"T + m' + (n''-np erfS + — e ^ (S'^'-S') 
1 c 1 Z /— 1 c. 

V TT 


continuity 


(U* +Up + <U^- U*) erfS + (n*- n|) ^ 

S J~r\ 
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... e be 
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1 2 s/S /S 1 2 s/S 1 2 


temperature 


-(S* +sp - (S*-Sp erfS + — (S*- s!) +— (S* - S* 
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-s 2 -s 2 -s 2 

( n T - np — (up Up — (n ‘ - np 

1 2 S/"n 1 2 SA^ 3/^ 1 2 


+ | ^ — (spsp + — (s*- s p = o 

3 1 2 3A7 12 


were used to simplify the above equations. Hence, in combining 
equations (D.ll) and (D. 13), the x-m omen turn equation becomes 


nu 2 (l + — 


1_ |u d 
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• 3 -\ 2 > 




(D. 16) 
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A similar expression can be obtained for the energy equation, 
but it is less obvious. Define 


X = 


o 


Y ■ s( 2+ V^4> 3 


2S 


Z = 4-^2+— + 


dx 


2 ' 4 

S 2S 


> 


The differential equations (D. 13) to (D. 15) simplify to 
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Au 
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where Q(S) = — . u fc 


-S 


. *'•' „ 'i' 




U 


- , AS = S - S_, AN = nT - nl - , and 

S/17 12 12 


Au - U'- - U . These three equations are solved for AS, AN, and Au 
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in terms of X, Y, and Z. The solution is 


A u 
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By placing these results into the energy equation (D. 12), we have 
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Equations (D. 16) and (D. 17) do not reduce to the one dimen- 
sional viscous equations even though they are similar in form. 
,22 

d/dx (2 + 3/S )u is related to a stress term and 
2 4 3 

d/dx (2 + 8/S + 5/2S )u can be thought of as a conduction type 
term. If the d/dx [ ] terms are zero, then (D. 10), (D. 16) and 

{D. 17) reduce to the inviscid equations 


nu = 
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or 

2 

p . u + P = Bi, 

.2 

PU ‘ (t - + 2 RT ) = B 3 

p u = Bj . 


The stress term in equations (D. 16) and (D. 17) is not invariant 
under a Galilean transformation as is the stress term in the one- 
dimensional viscous equations. This not surprising since Lees' two 
sided Maxwellian is referenced to a particular coordinate system. 



